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ABSTRACT 

Obtaining  far- field  patterns  in  electromagnetics  or  acoustics,  al¬ 
though  generally  not  as  computationally  expensive  as  solving  for 
the  sources  induced  on  an  object,  can  none-the-lessat  times  be  a 
substantial  fraction  of  the  overall  computer  time  required  for 
some  problems.  This  can  be  especially  the  case  in  determining 
the  monostatic  radar  cross  section  of  large  objects,  since  the  cur¬ 
rent  distribution  must  be  computed  for  each  incidence  angle,  or 
when  computing  the  radiation  patterns  of  large  reflector  antennas 
using  physical  optics.  In  addition,  when  employing  the  point 
sampling  and  linear  interpolation  of  the  far  field  that  is  most 
often  used  to  develop  such  patterns,  it  can  be  necessary  to  sam¬ 
ple  very  finely  in  angle  to  avoid  missing  fine  details  such  as 
nulls.  A  procedure  based  on  model-based  parameter  estimation 
is  described  here  that  offers  the  opportunity  of  reducing  the  num¬ 
ber  of  samples  needed  while  developing  an  easily  computed  and 
continuous  representation  of  the  pattern.  It  employs  windowed, 
low-order,  overlapping  fitting  models  whose  parameters  are  esti¬ 
mated  from  the  sparsely  sampled  far-field  values.  The  fitting 
models  themselves  employ  either  discrete-source  approximations 
to  the  radiating  currents  or  Fourier  models  of  the  far  field.  For 
the  cases  investigated,  as  few  as  1.5  to  2  samples  per  far-field 
lobe  are  found  to  be  sufficient  to  develop  a  radiation-pattern  esti¬ 
mate  that  is  accurate  to  0.1  dB,  and  2.5  samples  per  lobe  for  a 
simple  scatterer.  In  general,  however,  the  sampling  density  is 
not  determined  by  the  lobe  count  alone,  but  by  the  effective  rank 
of  the  field  over  the  observation  window,  which  in  turn  is  a 
function  of  both  the  aperture  size  and  the  spatial  variation  of  the 
source  distribution  within  that  aperture. 


1.  INTRODUCTION 

An  important  goal  of  all  numerical  modeling  is  that  of  mini¬ 
mizing  the  number  of  samples  needed  of  the  relevant  observables 
and  equations  so  as  to  minimize  the  total  computer  operation 
count  (OC)  (or  the  computer  cost)  while  achieving  a  desired  ac¬ 
curacy,  or  equivalently  reducing  the  uncertainty  to  a  specified 
level,  in  the  computed  results.  This  topic  is  considered  below  in 
the  context  of  computing  radiation  and  scattering  patterns  in 
electromagnetics.  The  goal  of  particular  interest  here  concerns 
minimizing  the  number  of  far-field  pattern  samples  that  are  re¬ 
quired  to  represent  a  radiation  pattern  and/or  the  number  of  inci¬ 
dence  angles  that  are  required  to  develop  a  monostatic  backscatter 
radar-cross  section  (RCS)  pattern.  An  added  benefit  of  the  proce¬ 
dure  described  below  is  that  of  obtaining  an  estimate  of  the  un¬ 
certainty  in  the  pattern  that  is  developed. 

The  approach  taken  employs  model-based  parameter  estimation 
(MBPE)  [1].  This  is  a  procedure  that  uses  reduced-order, physi¬ 
cally  based  fitting  models  (FMs)  whose  parameters  are  computed 
from  samples  of  first-principles  generating  models  (GMs)  such 
as  Maxwell  s  Equations.  Computation  of  a  FM  sample  is  trivial 
compared  with  one  evaluated  from  a  GM  for  large,  complex 
problems,  potentially  requiring  orders-of-magnitude  less  comput¬ 
ing  time.  This  makes  it  possible,  using  the  FM,  to  develop  an 


essentially  continuous  representation  of  physical  observables  of 
interest  as  opposed  to  the  pointwise  characterization  that  is  usu¬ 
ally  accepted  when  the  cost  of  eachGM  evaluation  is  large.  A 
list  of  acronyms  used  in  the  article  follows  the  references. 


2.  BACKGROUND 

One  of  the  most  frequently  encountered  problems  in  electromag¬ 
netic  field  computations  is  that  of  determining  a  radiation  or 
scattering  pattern  from  a  current  distribution(s)  known  over  some 
surface.  For  antenna  applications,  usually  a  single  current  distri¬ 
bution  is  of  interest,  while  for  RCS  computations,  a  new  current 
distribution  arises  for  each  incidence  angle  of  a  plane-wave  excit¬ 
ing  field.  In  either  case,  the  far  field  is  usually  needed  at  enough 
points  to  develop  a  smooth,  in  effect  continuous,  approximation 
of  the  overall  pattern  in  one  or  more  planes.  The  required  num¬ 
ber  of  radiation-pattern  samples  and  resulting  OC  are  proportion¬ 
al  to  the  maximum  body  dimension,  L,  in  the  plane  in  which 
the  pattern  is  being  computed.  Furthermore,  the  number  of  cur¬ 
rent  samples  on  the  surface,  S,  is  itself  is  proportional  to  the 

body  s  surface  area,  i.e.,  S  «  L2.  As  L,  and  therefore  S,  increas¬ 
es,  the  far-field  computation  can  become  a  significant  cost  in  ob¬ 
taining  radiation  and  scattering  patterns.  Thus,  reducing  the 
number  of  angle  samples  could  be  worthwhile  in  terms  of  reduc¬ 
ing  the  overall  computer  cost  of  obtaining  the  radiation  pattem(s) 
of  a  large  antenna  or  RCS  of  a  large  scatterer 

Some  previous  work  by  the  author  [2,3]  and  others  [4,5]  has  de¬ 
scribed  an  approach  that  uses  MBPE  to  decrease  the  number  of 
samples  that  are  needed  to  determine  a  far-field  pattern.  The 
work  presented  in  [2,3],  is  briefly  summarized  and  extended 
here,  especially  with  respect  to  how  the  errors  in  the  pattern  can 
be  estimated  and  controlled  and  how  the  pattern  itself  is  mod¬ 
eled.  Further  examples  of  efficient  pattern  computation  can  be 
found  in  the  work  of  Bucci  and  his  various  collaborators  who 
have  developed  signal-processing-like  techniques  for  computing 
far-field  patterns  [see  for  example  6,  7,  and  8]. 

3.  CHOOSING  THE  QUANTITY  TO  MODEL 

Estimating  a  far-field  pattern  using  MBPE  requires  choosing  the 
kind  of  reduced-order  FM  that  is  to  be  used  and  to  what  observ¬ 
able  that  FM  is  to  be  applied.  Two  obvious  choices  for  mini¬ 
mizing  the  number  of  pattern  samples  or  incidence  angles  are 
available:  1)  to  model  the  current  distribution;  or  2)  to  model  the 
far-field  pattern.  Furthermore,  there  are  two  ways  in  which  the 
current  distribution  itself  might  be  modeled.  These  various 
choices  are  briefly  summarized  below. 

3.1  Modeling  a  Current  (or  Aperture)  Distribution 

The  two  approaches  to  be  discussed  here  for  modeling 
the  current  (or  aperture)  can  be  best  described  as  Discrete-Source 
Approximations  (DSAs).  The  current  over  the  surface  of  the  ob¬ 
ject  under  consideration  or  field  over  an  aperture  is  replaced  by  a 
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linear  array  of  discrete,  or  point,  sources  aligned  in  space  as  is 
described  below.  (Other  point-source  geometries  might  also  be 
used;  attention  is  limited  here  to  a  linear  geometry).  The  param¬ 
eters  of  the  DSA  model  are  obtained  by  fitting  it  to  far-field 
samples  obtained  from  the  usual  integration  of  the  actual  current 
distribution.  The  DSA  is  then  used  to  approximate  the  pattern 
between  these  samples  to  thereby  obtain  a  continuous  estimate  of 
the  pattern  without  requiring  additional  current  integrations. 

3.1.1  Using  A  Prony  Model 

Prony  s  method  (or  its  equivalent)  can  be  used 
for  the  DSA  computation,  using  as  a  FM  [9] 

P 

F(0)  =  ^Saeikd«cos(e) 

a=l  (1) 

which  involves  P  point  sources  of  strengths  Sa  located  at  posi¬ 
tions  da  along  the  axis  of  the  DSA0  array,  with  0  the  angle  to 

the  far  field  measured  from  that  axis.  There  being  2N  unknown 
parameters  in  Eq.  1,  the  Prony  DSA  (PDSA)  thus  requires  2N 
far-field  samples  of  the  actual  pattern.  Furthermore,  these  2N 
samples  need  to  be  spaced  uniformly  in  cos(0),  a  feature  that 
makes  PDSA  less  suitable  for  the  adaptive-sampling  approach 
described  next.  A  possible  advantage,  however,  of  finding  the 
locations  of  the  discrete  sources  is  the  possibility  of  developing 
an  approximate  image  of  the  source  whose  far  field  is  being  sam¬ 
pled. 

There  are  two  different  candidate  DSA  geometries  that  might  be 
considered.  In  the  first,  shown  in  Fig.  la,  the  DSA  axis  is  fixed 
and  a  sequence  of  angle  windows  are  rotated  about  this  axis  over 
the  complete  range  of  observation  angles  of  interest.  In  the  sec¬ 
ond,  shown  in  Fig.  lb,  the  DSA  axis  is  itself  rotated  to  be  per¬ 
pendicular  to  the  angle  that  defines  the  midpoint  of  each  angle 
window  used  for  the  successive  DSA  computations.  Note  that 
L  varies  in  proportion  to  the  length  of  the  object  as  seen  from 
the  center  of  the  observation  window  when  using  the  approach  of 
Fig.  lb  whereas  L  is  fixed  at  the  maximum  linear  dimension  of 
the  object  in  the  observation  plane  for  approach  1  a. 

3.1.2  Using  A  Specified  DSA 

The  specified  DSA  (SDSA)  model  is  the  same 
as  Eq.  1  except  that,  since  the  source  locations  are  now  specified, 
only  the  N  source  strengths  Sa  are  unknowns.  As  for  the  Prony 

DSA,  samples  of  the  actual  far-field  are  used  to  obtain  the 
discrete-source  strengths.  In  contrast  to  the  Prony  model,  how¬ 
ever,  the  pattern  samples  are  not  constrained  in  their  placement 
but  can  be  arbitrarily  located  in  angle  and  only  N  are  needed,  two 
distinct  over  the  Prony  model.  For  the  SDSA  results  presented 
here,  the  sources  are  equally  spaced  along  the  array  axis,  with  the 
source  numbers  1  and  N  located  at  the  ends  of  the  aperture  L 
using  the  configuration  shown  in  Fig.  la.  As  the  order,  i.e., 
number  of  sources  used  in  a  particular  SDSA  FM,  is  increased, 
the  source  spacing  is  therefore  systematically  decreased  in  pro¬ 
portion  to  the  number  of  sources  that  are  used. 

3.2  Modeling  the  Pattern 

An  alternative  to  using  a  DSA  for  the  pattern  estimation 
is  to  model  the  pattern  itself,  using  a  Fourier  series,  an  approach 


denoted  as  the  Fourier  Series  Pattern  Model  (FSPM).  In  this 
case  the  FM  can  be  developed  as  [2] 

F(e)=^RaeM+  ^R'_ae-iae 

a  =  S  a^S*  (2) 

where  we  set  F  equal  to  Int(L  +  1),  with  Int(X)  denoting  the 
value  of  X  rounded  off  to  the  nearest  lower  integer.  Also,  2F  - 
S  -  S  =  N  with  N  the  total  number  of  terms  in  the  FM,  S  f  S 

f  S  +  1,  and  Ra  and  R  _a  are  the  amplitudes  of  the  positive  and 


Figure  1.  Two  ways  of  implementing  the  Discrete-Source  Array 
Fitting  Model.  In  (a)  the  sampling  window  (denoted  by  alternately 
light  and  heavy  lines)  rotates  in  angle  about  the  long  axis  of  the 
object  while  the  DSA  FM  axis  remains  fixed  along  which  the  dis¬ 
crete  sources  are  located.  In  (b),  on  the  other  hand,  the  sampling 
aperture  L’  is  rotated  with  respect  to  the  long  axis  of  the  object 
and  the  sampling  window  is  bisected  by  a  line  perpendicular  to  it. 
Forthe  case  of  the  PDSA,  an  N-source  FM  requires  a  minimum  of 
2N  field  samples.  Both  the  source  locations  (indicated  by  the 
X’s)  and  their  strengths  are  determined  by  sampling  the  far  field 
as  a  function  of  cos(0),  where  generally  the  sources  will  be  non- 
uniformly  spaced.  Alternatively,  for  the  SDSA  used  to  obtain  the 
results  presented  here,  the  configuration  (a)  was  used  with  N 
sources  uniformly  spaced  along  the  array  axis,  thus  requiring 
only  N  field  samples  for  computation  of  their  strengths. 
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negative  modes  respectively.  As  for  the  specified  DSA,  N  far- 
field  samples  are  required  to  quantify  the  parameters  of  the 
FSPM.  Note  that  the  Fourier  approach  yields  estimates  of  the 
far-field  Fourier  components  directly,  as  contrasted  with  either  of 
the  DSA  FMs  where  the  source  strengths,  and  locations  as  well 
for  the  PDSA,  are  the  parameters  being  computed.  Ultimately, 
however,  the  differences  between  the  SDSA  and  the  FSPM  are 
rather  slight,  differing  most  in  the  former  having  an  observation- 
angle  dependence  given  by  exp[ikdcos(0)]  while  the  latter  de¬ 
pends  on  exp(ia0). 

3.3  Adaptive  Estimation  Using  Windowed, 

Overlapping  Fitting  Models 

To  be  effectively  implemented,  any  adaptive  numerical 
procedure  requires  an  error  estimate.  For  the  specific  problem  of 
adaptive  pattern  estimation,  the  error  estimate  is  used  to  deter¬ 
mine  whether  a  given  FM  is  accurate  enough,  and  if  not,  where 
in  angle  a  new  sample  should  be  located.  Various  FM  configu¬ 
rations  might  be  considered  in  this  application.  For  example,  a 
single  FM  might  be  used  to  span  the  entire  angle  window  of  in¬ 
terest  with  an  initial  set  of  sample  angles,  S,  chosen  for  GM 
evaluation.  A  subset  of  these  GM  samples,  S-l,  could  be  used  to 
obtain  the  S-l  parameters  of  FM(S-l)  while  all  S  GM  samples 
could  be  similarly  used  to  find  the  S  parameters  of  FM(S). 
These  two  FMs  could  then  be  sampled  more  finely,  by  a  factor 
of  10  or  so,  in  angle  than  was  used  for  the  initial  GM  sampling, 
with  the  difference  between  them  serving  as  an  error  measure.  If 
the  error  measure  exceeds  the  allowable  uncertainty  specified  by 
the  modeler,  the  S+l  th  GM  sample  would  be  added  at  the  angle 
where  the  difference  between  FM(S-l)  and  FM(S)  is  greatest. 
The  parameters  of  FM(S+1)  could  then  be  computed  and  the  dif¬ 
ference  between  FM(S)  and  FM(S+1)  obtained.  The  process 
would  be  systematically  continued  until  the  maximum  difference 
between  FM(S+k-l)  and  FM(S+k)  satisfies  the  uncertainty  speci¬ 
fication,  with  S+k  the  total  number  of  GM  samples  required. 
The  initial  number  of  GM  samples  to  be  used,  S,  would  be  cho¬ 
sen  to  be  proportional  to  the  number  of  pattern  lobes  expected 
over  the  observation  window. 

Using  a  single  FM  to  cover  an  entire  pattern  would  generally  not 
be  the  best  approach,  however,  partly  due  to  the  growing  cost  of 
the  FM  computation  itself,  but  more  importantly  due  to  the  fact 
that  the  condition  number  of  the  data  matrix  needing  solution  for 
the  FM  parameters  may  increase  unacceptably,  especially  for  a 
large  number  of  pattern  lobes.  Instead,  as  is  done  here,  a  number 
of  windowed,  lower-order  (i.e.,  fewer  parameters)  overlapping 
FMs  are  used.  Each  FM  shares  two  or  more  GM  samples  with 
its  neighbors,  as  is  illustrated  conceptually  in  Fig.  2  where  there 
are  a  total  of  N  FMs.  As  above,  after  their  parameters  have  been 
evaluated,  the  FMs  are  evaluated  more  finely  in  angle  than  was 
the  GM  pattern  initially,  and  the  differences  between  the  sets  of 
overlapping  FMs  are  computed.  A  new  GM  sample  then  is 
added  where  the  maximum  difference  between  all  sets  of  overlap¬ 
ping  FMs  is  found  to  occur.  The  respective  parameters  of  the  af¬ 
fected  FMs  (i.e.,  those  whose  windows  contain  the  new  sample 
angle)  are  then  updated,  with  the  process  of  computing  FM  dif¬ 
ferences  and  newGM  samples  continuing  until  the  specified  un¬ 
certainty  is  satisfied  over  the  entire  pattern  window. 

An  appropriate  choice  for  the  number  of  FMs  for  a  given  prob¬ 
lem  might  require  some  experimentation.  If  it  s  desired  to  keep 
the  order  of  all  FMs  below  some  specified  value,  then  based  on 
the  computer  experiments  done  in  getting  the  results  presented 
below,  eachFM  should  have  an  angle  window  that  spans  just  a 


few  pattern  lobes,  say  three  or  four  at  most.  The  number  of  lobes 
to  be  expected  can  generally  be  estimated  from  the  size  of  the 
aperture  whose  pattern  is  being  modeled.  However,  it  may  hap¬ 
pen  that  adding  a  new  GM  sample  to  a  FM  causes  its  to  rank  ex¬ 
ceed  some  specified  limit.  By  dividing  such  a  FM  into  two 
lower-order,  overlapping  ones,  the  problem  of  excess  FM  rank 
can  be  avoided  and  the  initial  number  of  FMs  will  not  be  so  im¬ 
portant. 


Figure  2.  A  conceptual  illustration  of  the  use  of  overlapping  fitting 
models  for  adaptive  sampling  of  a  far-field  pattern.  The  horizontal 
lines  indicate  the  angular  extent  of  each  FM  with  the  open  circles 
showing  where  the  pattern  (or  GM)  is  initially  sampled.  In  this  par¬ 
ticular  case,  all  interior  FMs  begin  with  three  GM  samples  while 
those  at  either  end  use  only  two,  for  a  total  of  N  + 1  starting  sam¬ 
ples.  Additional  GM  samples  are  systematically  located  where 
the  maximum  difference  occurs  between  two  overlapping  FMs 
until  a  specified  convergence  criterion  has  been  satisfied.  Only 
two  FMs  overlap  in  any  of  the  observation  windows  here,  but 
other  overlap  arrangements  could  be  used  as  well. 

3.4  Specification  of  the  Fitting-Model  Error  and 

Computing  the  Final  Pattern  Estimate 

In  contrast  to  specifying  a  fixed  FM  difference,  or 
fitting-model  error  (FME),  between  overlapping  FMs  as  was 
done  in  [2,3],  the  maximum  FME  can  also  be  scaled  relative  to 
how  the  magnitudes  of  the  far-field  samples  vary.  For  example, 
as  the  magnitude  decreases  the  FME  might  be  proportionately 
increased  to  accommodate  the  fact  that  side-lobe  maxima  may 
not  be  needed  to  the  same  accuracy  as  the  main  lobe.  In  the  re¬ 
sults  to  follow,  the  maximum  permitted  FME  at  the  observation 
angle  0  is  given  by 


a  =  1 

(3) 

where  Aj  and  A2  are  specified  parameters,  GM(0max)  is  the 
maximum  value  of  the  pattern  being  modeled,  and  FMa(0)  is 

the  value  of  the  a  th  FM  at  0  where  a  total  of  M  FMs  overlap. 
The  parameter  Aj  determines  the  maximum  acceptable  FME  in 
the  vicinity  of  the  peak(s)of  the  GM.  The  parameter  A2  increas¬ 
es  the  allowable  FME  in  proportion  to  the  decrease  in  the  FM 
values  relative  to  GM(0max).  For  the  results  that  follow,  the 

nominal  values  for  Aj  and  A2  were  0.1  and  0.05  dB,  respective- 


FME(0)  =  Ai  + 
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ly ,  but  FME(0)  limited  to  a  maximum  value  of  3  dB.  There  are 
many  possible  variations  that  might  be  used  for  setting  the  ac¬ 
ceptable  error  or  uncertainty,  among  which  would  be  relaxing  the 
desired  accuracy  in  certain  angular  regions  or  increasing  it  in  oth¬ 
ers  consistent  with  the  requirements  of  a  particular  application. 

The  final  pattern  estimate  is  then  obtained  by  computing  the  av¬ 
erage  of  the  successively  overlapping  FMs  as  the  observation 
angle  is  scanned  over  the  angle  pattern.  Thus,  for  M  overlapping 
FMs  at  angle  0  we  would  have 

Fave,M(0)  =  j^[Fi(6)  +  Fi+l(9)  "*  +  Fi+M-ll. 

(4) 

For  the  results  presented  below,  M  =  2  was  used. 

Note  that  in  addition  to  controlling  the  adaptation  process,  the 
error  specification  provides  a  measure  of  the  uncertainty  in  the 
average  FM  values  of  Eq.  4.  Since  Eq.  3  gives  the  maximum 
variation  permitted  between  the  FMs  that  overlap  at  a  given 
angle,  it  s  proposed  that 


window  to  be  half  overlapped  with  its  nearest  neighbors 
seems  to  be  a  good  choice.  Note  that  using  half- 
overlapped  FMs  leads  naturally  to  a  minimum  of  3  ini¬ 
tial  GM  samples  per  FM  window  as  shown  in  Fig.  2. 

6)  Finally,  the  kind  of  FM  error  to  be  used  and  its  nu¬ 
merical  parameters,  as  in  Eq.  (3),  must  be  specified. 


4.  USING  THE  SPECIFIED  DISCRETE-SOURCE 
APPROXIMATION 

In  previous  work  [2,3]  the  FSPM  approach  was  described  and 
some  initial  results  presented,  demonstrating  that  the  far-held 
pattern  could  be  estimated  to  an  amplitude  uncertainty  of  0.01 
(about  0.05  dB)  using  only  ~  3  to  3.5  samples  per  lobe.  Several 
results  obtained  using  the  alternate  SDSA  model  are  presented 
here  together  with  an  illustration  of  how  the  model  performance 
depends  on  the  uncertainty  specification.  The  particular  patterns 
used  for  testing  the  MBPE  approach  here  were  chosen  as  being 
representative  of  the  kinds  of  patterns  encountered  in  typical  ap¬ 
plications,  as  well  as  their  having  closed-form  analytical  expres¬ 
sions,  except  for  the  last  example  of  the  random-source  array. 


Fave,M±(0)  =  Fave,M(0)  ±  ”FME(0) 

1  (5) 

will  yield  realistic  upper  and  lower  error-bound  estimates  for  the 
average  pattern  values.  That  Eq.  5  does  indeed  provide  a  realis¬ 
tic  error  bound  relative  to  the  FM  estimates  of  GM  is  demon¬ 
strated  in  the  following  examples. 

It  s  worth  emphasizing  that,  in  contrast  to  using  an  entire- 
domain  basis  for  the  far  field,  cylindrical-  or  spherical-wave  ex¬ 
pansions,  for  example,  the  windowed  approach  employed  here 
together  with  discrete  point  sources  can  be  applied  to  smaller,  or 
limited,  angular  sectors  with  no  computational  penalty.  Also 
observe  that  in  implementing  the  SDSA,  angle  sampling  is  done 
in  terms  of  cos(0)  rather  than  0.  This  is  because  the  pattern 
lobes  tend  to  be  distributed  uniformly  in  terms  of  the  former 
variable. 

To  summarize  then,  the  following  steps  are  involved  in  the 
MBPE  modeling  of  radiation  and  scattering  patterns: 

1)  The  quantity  to  be  modeled,  i.e.,  the  source  distribu¬ 
tion  or  pattern  itself,  is  first  selected. 

2)  The  total  angle  window  over  which  the  pattern  is  to 
be  estimated  is  specified. 

3)  The  number  of  pattern  lobes,  L,  that  are  anticipated 
over  the  angle  window  of  interest,  considering  the  aper¬ 
ture  size  in  wavelengths  whose  pattern  is  to  be  found,  is 
estimated. 

4)  The  number,  N,  of  initial  fitting  models  to  be  used 
then  needs  to  be  chosen.  Choosing  a  value  for  N  be¬ 
tween  L/2  and  2L/3,  with  L  expressed  in  wavelengths, 
should  provide  a  reasonable  starting  point,  noting  that 
the  smaller  the  value  of  N  that  is  used  the  larger  will  be 
the  required  FM  order. 

5)  The  configuration  of  the  FM  overlap  then  needs  to 
be  selected.  This  can  be  fairly  flexible.  Arranging  each 


4.1  The  Far  Radiated  Field  of  the  Uniform  Current 
Filament 

The  far  field,  F(0)ucp,  of  a  uniform  current  filament  (or 
uniform  aperture)  is  proportional  to  [10] 


F(0)ucf  —  L 


sin[7iLsin(0)] 

7rLsin(0) 


(6) 


where  the  filament  length  L  is  expressed  in  wavelengths  and  the 
observation  angle  0  is  measured  from  the  filament  axis.  A  total 
of  N  =  13  FMs  was  used  for  an  L  =  20-wavelength  UCF,  ar¬ 
ranged  as  shown  in  Fig,  2,  with  the  parameters  of  Eq.  3  given 
by  Ay  =0.1  and  A2  =  0.05  dB.  Thus,  all  FMs  initially  have 
three  GM  samples  except  for  those  on  either  end,  which  have 
only  two.  This  results  in  eachFM  sharing  two  samples  with  its 
nearest  neighbors  except  for  those  on  the  end,  which  overlap 
with  one  adjacent  FM.  Applying  the  SDSA  to  GM  samples  of 
Eq.  6  yields  the  results  of  Fig.  3,  normalized  to  a  maximum  of 
0  dB,  where  the  upper  and  lower  error-bound  estimates  for  the 
pattern  peaks,  using  Eq.  5,  and  the  actual  pattern  from  Eq.  6  are 
plotted.  The  actual  pattern  is  seen  to  lie  between  the  upper-  and 
lower-bound  peaks  throughout  the  entire  window.  In  Fig.  4  the 
average  values  of  the  overlapping  FMs  from  Eq.  4,  FMave,  are 

compared  with  the  actual  pattern,  where  34  of  the  35  GM  sam¬ 
ples  used  for  computing  the  final  FM  parameters  are  also  shown. 
The  maximum  difference  between  the  actual  pattern  and  FMave 

is  consistent  with  the  error  specification  of  Eq.  4  and  the  numeri¬ 
cal  values  used  for  its  parameters.  With  35  samples  needed  to 
estimate  a  pattern  having  20  lobes  or  maxima,  1.75  samples  are 
required  per  lobe  for  the  UCF. 

The  behavior  of  three  different  error  measures  for  the  L  =  20 
UCF  is  presented  in  Fig.  5  as  a  function  of  the  number  of  GM 
samples  used  for  computing  the  FM  parameters.  The  upper  plot 
on  this  graph  shows  the  maximum  difference  between  all  pairs  of 
the  13  finely  sampled,  overlapping  FMs  as  a  function  of  the 
number  of  GM  samples  used  for  their  computation.  The  middle 
plot  exhibits  the  difference  between  all  pairs  of  overlapping  FMs 
averaged  over  the  angles  they  commonly  sample.  The  bottom 
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plot  displays  the  angle-averaged  difference  between  FMave  and 

the  actual  pattern,  or  the  GM.  Of  these  three,  only  the  first  two 
would  be  available  in  actual  application  since  the  GM  samples 
needed  to  compute  the  last  error  measure  would  be  limited  to 
those  available  up  to  that  point  in  the  adaptation  process.  It  s 
useful  to  include  the  latter  error  measure,  however,  as  the  differ¬ 
ence  between  the  FMaye  and  the  GM  provides  a  reality  check  on 

the  MBPE  performance. 


Observe  that  the  difference  between  the  overlapping  FMs  for  this 
example  always  exceeds  the  corresponding  FMave  -  GM  differ¬ 
ence.  This  shows  that  the  FM-FM  error  measure  provides  a  con¬ 
servative,  or  high-side,  estimate  of  the  error  (or  uncertainty)  in 
FM„VP  relative  to  the  true  GM  values.  Note  also  that  the  maxi- 

mum  FM  difference  does  not  decline  monotonically.  This  is  be¬ 
cause  updating  two  FMs  can  at  times  increase  the  maximum  dif¬ 
ferences  that  then  results  with  their  overlapping  neighbors. 


COSINE  OF  OBSERVATION  ANGLE 


Figure  3.  One  quadrant  of  the  normalized 
radiation  pattern  for  the  field  of  an  L  =  20 
(wavelengths)  UCF  as  obtained  using  the 
specified  discrete-source  approximation. 
The  lower- and  upper-bound  peaks  estimat¬ 
ed  from  Eq.  5  are  shown  together  with  the 
actual  pattern  from  Eq.  6.  The  pattern  is  cut¬ 
off  at  -40  dB  to  show  the  difference  be¬ 
tween  the  various  plots  more  clearly. 
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Figure  4.  One  quadrant  of  the  normalized,  D 
actual  field  of  the  20-wavelength  UCF  com-  ^ 
pared  with  FMave.  The  circles  indicate  34  "20 

of  the  35  GM  samples  computed  from  Eq.  6  q; 

(one  is  below  -60  dB),  that  were  used  in  the  <  -30 
estimation  process.  Because  the  cosine-  u" 
angle  values  used  to  obtain  the  average  Q 
and  actual  patterns  do  not  always  precisely  "40 
match  those  used  to  obtain  the  original  GM  -s 
samples,  some  of  the  latter  do  not  coincide  <  .gQ 
with  the  pattern  plots.  S 
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NUMBER  OF  GM  SAMPLES 


Figure  5.  Three  error  measures  discussed 
in  the  text  for  the  UCF  shown  as  a  function 
of  the  number  of  GM  samples  used  up  to 
that  point  in  the  adaptation  process  over  an 
observation  window  *90  to  0  deg,  with  the 
parameters  in  Eq.  3  given  by  A-)  =  0.1  dB 

and  A2  =  0.05  dB.  The  open  squares  de¬ 
note  the  maximum  pointwise  difference  be¬ 
tween  all  pairs  of  overlapping  FMs  with  the 
addition  of  a  new  GM  sample.  The  closed 
and  open  circles  exhibit  the  angle-averaged 
difference  between  the  overlapping  FMs 
and  between  the  average  FM  and  the  GM  it¬ 
self.  The  fact  that  the  latter  is  smaller  than 
the  former  shows  that  the  FM  difference 
provides  a  conservative  measure  of  the  ac¬ 
curacy  of  FMave. 


Figure  6.  One  quadrant  of  the  normalized 
radiation  pattern  for  the  field  of  an  L  =  20 
(wavelengths)  SCF  as  obtained  using  the 
SDSA.  The  lower- and  upper-bound  peaks 
estimated  from  Eq.  5  are  shown  together 
with  the  actual  field  from  Eq.  7.  The  pattern 
is  plotted  over  a  limited  dB  range  to  show 
the  difference  between  the  various  plots 
more  clearly. 


COSINE  OF  OBSERVATION  ANGLE 


COSINE  OF  OBSERVATION  ANGLE 


Figure  7.  One  quadrant  of  the  actual  field  of 
the  20-wavelength  SCF  compared  with 
FMave.  The  circles  indicate  28  of  the  40  GM 

samples,  computed  from  Eq.  7,  that  were 
used  in  the  estimation  process. 
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4.2  The  Far  Radiated  Field  of  the  Sinusoidal 
Current  Filament 

The  normalized  far  field  of  a  center-fed  sinusoidal  cur¬ 
rent  filament  (SCF)°L  wavelengths  long  is  proportional  to  [10] 


Fscf(G) _ 


cos(7iLcos9)  -  cos(ttL) 
sin0 


(7) 


Application  of  adaptive  MBPE  to  an  L  =  20-wavelength  SCF 
using  N  =  13  overlapping  FMs  as  shown  in  Fig.  2  leads  to  the 
upper-  and  lower-bound  estimates  for  the  far  field  shown  together 
with  the  actual  pattern  in  Fig.  6.  Again,  the  actual  field  values 
are  seen  to  lie  between  the  bounding  curves  provided  by  MBPE 
adaptive  sampling.  A  comparison  of  FMave  with  the  actual  pat¬ 
tern  is  presented  in  Fig.  7  where  most  of  the  40  GM  samples 
used  for  the  FM  computation  are  also  shown  (some  fall  below  - 
20  dB).  The  actual  and  FMave  curves  are  essentially  graphically 

indistinguishable  on  the  scale  used.  The  20-wavelength  SCF 
has  only  10  lobes  ratherthan  the  20  lobes  of  the  20-wavelength 
UCF  over  the  same  -1  to  0  cos(0)  interval,  but  requires  6  more 
unknowns  to  achieve  the  same  estimation  uncertainty,  or  4  sam¬ 
ples  per  lobe.  This  indicates  that  the  number  of  GM  samples 
needed  to  achieve  a  given  pattern-estimation  uncertainty  is  sensi¬ 
tive  not  only  to  the  lobe  structure  of  the  pattern  itself,  but  also 
to  the  spatial  variation  of  the  source  distribution  that  produces 
that  pattern. 


4.3  The  Far  Field  Scattered  from  a  Thin,  Circular 
Cylinder 

The  GM  here  is  the  approximate  far,  scattered  field  of  a 
thin,  circular  cylinder  L  wavelengths  long,  which  is  proportional 
to  [11] 


RCS~  Fcyl(0)  = 


sin^TiLsinBi) 

cosQi— - — — 

zTtLsinBi 


2 


(8) 


for  which  MBPE  adaptive  sampling  produces  the  results  shown 
in  Figs.  8  and  9  for  L  =  10  wavelengths  with  7  FMs  being  used. 
In  contrast  with  an  L-wavelength  SCF,  where  the  number  of 
lobes  over  a  -1  to  0  interval  in  cos(0)  is  ~  L/2  and  the  UCF 
which  has  ~  L  lobes  over  that  same  interval,  the  L-wavelength 
scattererhas  ~  2L  lobes,  20  in  this  example.  As  for  the  previous 
cases,  the  FMave,  as  obtained  from  48  GM  samples,  and  a  finely 

sampled  GM  plot  for  the  scatterer  are  graphically  indistinguish¬ 
able.  The  number  of  samples  required  per  lobe  for  the  10- 
wavelength  cylinder  scatterer  is  thus  2.4,  greater  than  the  1.75 
needed  for  the  UCF  even  though  their  patterns  are  quite  similar, 
as  can  be  seen  by  comparing  Figs.  4  and  9. 


4  4  par  Radiated  Field  of  a  Randomized  Array 
of  Point  Sources 

The  last  example  considered  here  is  a  linear  array  of  2 1 , 
isotropic  point  sources  having  random  amplitudes  varying  be¬ 
tween -1  and  +1  and  located  randomly  within  a  10-wavelength 
aperture.  Results  for  cos(0)  varying  from  -1  to  0  are  shown  in 
Figs.  10  and  11,  again  using  13  FMs.  A  total  of  26  GM  sam¬ 
ples  is  needed  to  achieve  the  same  specified  estimation  error  as 
used  for  the  previous  cases.  Since  there  are  on  the  order  of  8.5 
maxima  in  the  pattern,  this  works  out  to  about  3  samples  per 


lobe,  midway  between  that  required  for  the  scatterer  and  the 
SCF.  The  number  of  samples  needed  per  lobe  or  per  wavelength 
of  aperture  to  achieve  the  same  specified  estimation  accuracy  (as 
given  by  Eq.  3  with  A^  and  A2  0.1  dB  and  0.05  dB,  respective¬ 
ly)  for  the  various  sources  considered  here  is  summarized  in 
Table  I  below. 


5.  CONCLUDING  OBSERVATIONS  CONCERNING 
ESTIMATION  UNCERTAINTY  AND  SAMPLING 

The  preceding  examples  demonstrate  that  adaptive  sampling  of 
radiation  and  scattering  patterns  using  MBPE  with  discrete- 
source-approximation  FMs  can  be  effective  in  not  only  reducing 
the  number  of  samples  needed  to  obtain  a  reduced-order,  continu¬ 
ous  representation  of  a  pattern,  but  also  in  constraining  the  esti¬ 
mated  pattern  to  satisfy  an  uncertainty  specification.  Some  addi¬ 
tional  computations  are  included  here  to  shed  further  light  on  the 
sampling  requirements.  The  L  =  20  UCF  problem  was  repeated 
using  Aj  =  0.05  and  A2  =  0  dB  to  reduce  the  estimation  uncer¬ 
tainty  to  a  smaller,  and  constant,  value  compared  with  the  criteri¬ 
on  used  in  obtaining  the  previous  results.  Using  these  coeffi¬ 
cient  values  in  the  FME  given  by  Eq.  3  results  in  a  maximum 
acceptable  difference  between  overlapping  FMs,  and  hence,  a 
maximum  estimation  error,  of  no  more  than  0.05  dB.  This  is 
much  less  than  might  normally  be  sought  in  practice  but  pro¬ 
vides  a  more  stringent  test  of  the  MBPE  procedure.  It  should  be 
noted  that  if  the  GM  samples  are  of  limited  accuracy,  for  exam¬ 
ple  being  derived  from  a  numerical  first-principles  model  for  a 
complex  problem,  then  seeking  an  accuracy  in  the  pattern  esti¬ 
mate  that  the  GM  samples  cannot  provide  might  result  in  stag¬ 
nating  the  estimation  process,  i.e.,  convergence  may  not  occur. 
But  when  using  an  analytical  expression  for  a  pattern,  as  is  done 
here,  this  will  not  be  a  problem.  It  s  also  worth  noting  that  the 
variational  relationship  between  the  far  fields  and  the  sources  that 
produce  them  results  in  errors  in  the  latter  not  translating  into 
comparable  errors  in  the  former. 

Upon  running  the  L=  20-wavelength  UCF  problem  using  these 
new  values  for  the  FME  coefficients  and  plotting  the  same  error 
measures  as  for  Fig.  5,  the  results  shown  in  Fig.  12  are  ob¬ 
tained.  Also  included  in  Fig.  12  are  lines  of  the  form  Aexp(-Bx) 
where  x  is  the  number  of  GM  samples  and  A  and  B  are  best-fit 
parameters.  It  s  interesting  to  see  that  the  three  best-fit  lines  are 
nearly  parallel,  with  all  decreasing  exponentially  as  a  function  of 
the  number  of  GM  samples.  This  behavior  is  similar  to  the  con¬ 
vergence  of  numerical  solutions  for  various  wire  geometries  as 
the  number  of  unknowns  in  a  moment-method  solution  is  in¬ 
creased.  A  0. 1  dB  in  the  angle-averaged  FM  difference  and  the 
difference  between  the  FMave  and  the  actual  pattern  is  achieved 
using  about  37  GM  samples. 

This  computation  was  repeated  for  UCFs  having  lengths  of  L  = 
10,  15  and  25  wavelengths,  yielding  the  results  of  in  Fig.  13 
where,  for  clarity,  only  the  FM-FM  differences  are  plotted.  The 
data  for  each  of  these  UCF°  lengths  are  found  to  be  fit  compara¬ 
bly  well  by  exponentials  having  different  slopes.  This  is  an  ex¬ 
pected  result  since  the  number  of  GM  samples  needed  to  achieve 
a  given  uncertainty  for  a  specified  source  distribution  is  expected 
to  be  related  to  the  number  of  pattern  lobes  which  are,  in  turn, 
proportional  to  the  aperture  size.  When  the  latter  effect  is  re¬ 
moved  by  replotting  the  data  of  Fig.  13  as  a  function  of  the 
number  of  GM  samples  per  wavelength  of  aperture,  the  results 
shown  in  Fig.  14  are  obtained.  The  best-fit  lines  for  the  various 
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UCF  lengths  are  nearly  coincident,  shifting  to  the  left  slightly 
with  increasing  L  and  exhibiting  slopes  that  are  within  about 
1.5%.  A  0.1  dB  angle-averaged  FM  difference  is  seen  to  require 
about  2  samples  per  wavelength  of  aperture. 


Since  a  UCF  L-wavelengths  in  extent  produces  L  lobes  per  90 
deg  in  its  pattern,  this  sampling  density  translates  to  about  2 
samples  per  lobe  as  well.  As  can  be  deduced  from  Table  I,  it  s 
clear  that  the  far-field  sampling  density  depends  on  more  than 
just  the  lobe  count  in  the  pattern. 
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Figure  8.  One  quadrant  of  the  normalized 
scattering  pattern  for  the  field  of  an  L  =  1 0 
wavelengths  thin  cylinder  as  obtained  using 
the  SDSA.  The  lower-  and  upper  -bound 
peak  estimates  of  the  far  field  from  Eq.  5  are 
shown  together  with  the  actual  field  from 
Eq.  8. 


Figure  9.  One  quadrant  of  the  actual  scat¬ 
tered  field  of  the  10-wavelength  cylinder 
compared  with  FMave.  The  circles  indicate 

some  of  the  48  GM  samples,  computed 
from  Eq.  8,  that  were  used  in  the  estimation 
process. 
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How  this  sampling  density  might  generalize  to  2D  or  3D  source 
distributions  needs  to  be  considered.  First  note  that  the  sam¬ 
pling  density  of  2  per  wavelength  of  aperture  for  the  UCF  ap¬ 
plies  to  a  pattern  symmetric  about  broadside  where  an  angle- 
observation  window  of  only  90  deg  needs  to  be  sampled.  For  a 
non-symmetric,  but  otherwise  equivalent,  linear  source  distribu¬ 
tion  this  would  then  work  out  to  4  samples  per  wavelength  since 
a  far- field  window  of  180  deg  would  then  need  to  be  sampled. 
For  a  true  2D  source  whose  pattern  has  to  be  sampled  over  360 
deg,  this  would  imply  that  8  samples  per  wavelength  of  its  max¬ 
imum  linear  aperture  are  needed  to  achieve  a  comparable  pattern 
accuracy  of  0.1  dB.  This  measure  would  seem  to  hold  in  planar 


cuts  for  the  radiation  patterns  of  3D  source  distributions  as  well. 

Results  obtained  for  the  cylinder  scatterer  indicate,  on  the  other 
hand,  that  nearly  5  samples  per  wavelength  of  cylinder  length  are 
required  for  the  same  kind  of  pattern-estimation  accuracy  for  a 
90-deg  sector  of  the  scattered  field  to  be  estimated.  Following 
the  same  line  of  reasoning  as  above  for  a  general  2D  scatterer,  or 
for  planar  cuts  of  3D  objects,  sampling  over  360  deg  can  be  esti¬ 
mated  to  require  on  the  order  of  20  samples  per  wavelength. 
One  obvious  question  arises  about  why  a  greater  sampling  densi¬ 
ty  is  apparently  needed  for  a  scattering  pattern  as  compared  with 
a  radiation  pattern  when  the  two  are  similar-appearing,  as  exhib- 
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ited  by  Figs.  4  and  9?  Perhaps  the  answer  is  that  for  the  latter 
situation  a  single  current  distribution  produces  the  entire  radiated 
field.  For  the  former,  on  the  other  hand,  the  current  distribution, 
being  a  function  of  the  angle  of  incidence,  changes  with  every 
viewing  angle  when  determining  the  monostatic  radar  cross  sec¬ 
tion. 

It  would  also  seem  to  follow  then  that  the  number  of  GM  sam¬ 
ples  needed  to  achieve  a  specified  estimation  uncertainty  must 
depend  on  both  the  aperture  size  and  spatial  variation  of  the 
source  within  it,  a  dependence  that  ultimately  is  exhibited  by  the 
pattern  rank  (R)  or  the  number  of  degrees  of  freedom  of  the 
pattern,  over  the  chosen  observation  window.  As  a  specific  ex¬ 
ample  of  how  the  source  characteristics  and  aperture  might  influ¬ 
ence  the  rank,  consider  an  aperture  L-wavelengths  long  having  N 
point  sources  uniformly  spaced  within  it. 


For  the  simplest  case  of  two  point  sources,  as  L  is  increased,  the 
number  of  pattern  lobes  over  360  deg  in  an  observation  plane 
containing  the  sources  will  be  of  order  4L.  However,  R,  as  de¬ 
termined  by  eigenvalue  analysis  of  the  data  matrix  that  arises 
when  using  Prony  s  method  remains  fixed  at  two  when  L  ex¬ 
ceeds  0.5,  regardless  of  how  large  the  aperture  is  made  and  how 
many  lobes  are  in  the  pattern.  On  the  other  hand,  if  L  is  fixed 
and  the  number  of  point  sources  is  systematically  increased,  R 
will  increase  proportionately  until  N  >  2L  whereupon  R  remains 
fixed  at  ~  2L  since  only  about  2  source/wavelength  can  be  re¬ 
solved  in  the  far  field.  Because  it  can  model  such  point-source 
arrays,  a  Prony-based  procedure  can  exploit  the  reduced  rank  of 
such  special  problems,  but  the  discrete-source  approximation  and 
a  Fourier  model  of  the  far  field  are  not  as  well-suited  for  doing 
so.  Generally  speaking,  with  everything  else  being  equal,  the 
best  pattern  estimator  would  be  one  for  which  the  number  of 
GM  samples  can  be  reduced  to  as  close  to  R  as  possible.  More 
investigation  is  needed  to  settle  these  issues,  and  to  generalize 
the  results  beyond  the  simple  cases  considered  here. 


Figure  10.  One  quadrant  of  the  radiation 
pattern  for  the  field  of  21  isotropic  point 
sources  having  random  amplitudes  be¬ 
tween  -1  and  +  1  and  randomly  located 
along  a  an  L  =  10  wavelengths  linear  array. 
Upper-bound  estimates  of  the  pattern  are 
shown  by  the  x’s  and  the  lower-bound  esti¬ 
mates  by  the  open  circles  with  the  actual 
pattern  shown  by  the  solid  line. 
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Figure  11.  The  actual  pattern  of  21  isotrop-  w 
ic  sources  having  random  amplitudes  and  O 
located  at  random  positions  along  a  10-  m  0 
wavelength  aperture  compared  with  £ 
FMave.  The  circles  indicate  the  26  GM  sam- 

pies  used  in  the  estimation  process.  < 
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Figure  12.  The  three  error  measures  shown 
in  Figure  6  are  plotted  here  for  the  20- 
wavelength  UCF  but  with  Ai  =  0.05  and  A2 

=  0  in  Eq.  3  to  obtain  an  estimate  that  is  both 
more  accurate  and  remains  constant  over 
the  dynamic  range  covered  by  the  pattern 
sampled  in  a  -90  to  0  deg  observation  win¬ 
dow.  A  best-fit  straight  line  is  computed  for 
all  three,  where  RA2  is  the  square  of  the  cor¬ 
relation  coefficient.  All  three  best-fit  lines 
are  seen  to  be  nearly  parallel  on  this  log- 
linear  plot,  being  well-described  as  expo¬ 
nential  functions  of  the  number  of  GM  sam¬ 
ples. 


100 

Figure  13.  The  angle-averaged  FM-FM  dif-  § 
ferences  for  a  UCF  varying  from  10  to  25  w  10 
wavelengths  using  A-|  =  0.05  and  A2  =  0  in  HI 

Eq.  3  yields  the  results  shown  here,  again  2 
for  a  -90  to  0  deg  observation  window.  HI  ' 
Best-fit  straight  lines  are  also  shown  togeth-  Jjj 
er  with  their  respective  correlation  coeffi-  u. 
cients.  The  constants  in  the  exponentials  It  -1 
for  the  best-fit  lines  are  seen  to  decrease  O 
with  increasing  filament  length.  The  expo-  *g 
nential  decrease  of  the  error  measures  as  U.  .01 
the  number  of  GM  samples  increases  is  a  g 
very  beneficial  attribute  of  a  numerical  pro-  LL 
cess.  .001 
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Figure  14.  Upon  replotting  the  results  of 
Fig.  13  but  with  the  number  of  GM 
samples/wavelength  as  the  independent 
variable,  the  best-fit  angle-averaged  FM-FM 
differences  become  nearly  coincident. 
These  results  show  that  for  this  particular 
case  the  number  of  GM  samples  needed  for 
a  specified  angle-averaged  FM  difference  in 
the  pattern  is  linearly  proportional  to  the  fila¬ 
ment  length.  They  also  show  that  the  sam¬ 
pling  density  required  for  a  given  value  of 
error  measure  decreases  somewhat  as  the 
UCF  increases  in  length. 
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TABLE  I:  FAR-FIELD  SAMPLING  REQUIRED  OVER  A  90- 
deg  WINDOW  FOR  SOURCES  TESTED  USING  Al  =  0.1 
AND  A2  =  0.05  dB  in  Eq.  3. 


SOURCE  TYPE 

NUMBER  FAR-FIELD 

FAR-FIELD 

TESTED 

OF  SAM- 

SAMPLES/ 

SAMPLES/ 

PLES 

LOBE 

WAVELENGTH 

OF  APERTURE 

20-WL 

SINUSOIDAL 

CURRENT 

FILAMENT 

40 

4 

2 

10- WL 

RANDOMIZED 

ARRAY 

26 

3 

2.6 

10- WL  CYLINDER 
SCATTERER 

48 

2.4 

4.8 

20-WL  UNIFORM 
CURRENT 
FILAMENT 

35 

1.75 

1.75 
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ABSTRACT:  In  this  paper  we  illustrate  how  to  solve 
the  general  Helmholtz  equation  starting  from 
Laplace's  equation.  The  interesting  point  is  that  the 
Helmholtz  equation  has  a  frequency  term  whereas 
Laplace’s  equation  is  the  static  solution  of  the  same 
boundary  value  problem.  In  this  new  formulation  the 
frequency  dependence  is  manifested  in  the  form  of  an 
excitation.  A  new  boundary  integral  method  for 
solving  the  general  Helmholtz  equation  is  developed. 
This  new  formulation  is  developed  for  the  two- 
dimensional  Helmholtz  equation  with  the  method  of 
moments  Laplacian  solution.  The  main  feature  of  this 
new  formulation  is  that  the  boundary  conditions  are 
satisfied  independent  of  the  region  node 
discretizations.  The  numerical  solution  of  the  present 
method  is  compared  with  finite  difference  and  finite 
element  solutions  of  the  same  problem.  Application  of 
this  method  is  also  presented  for  the  computation  of 
cut-off  frequencies  for  some  canonical  waveguide 
structures. 

L  INTRODUCTION 

The  two  dimensional  Helmholtz’s  equation  is  an 
important  equation  to  be  solved  in  many  numerical 
electromagnetics  problems  such  as  waveguide  related 
problems  and  appears  in  a  variety  of  physical 
phenomena  and  engineering  applications,  such  as, 
acoustic  radiation  [1],  heat  conduction  [2],  and  water 
wave  propagation  [3],  In  semiconductor  device 
modeling,  Helmholtz’s  equation  arises  frequently  as 
an  intermediate  step  in  the  solution  of  the  nonlinear 
Poisson’s  problem.  To  solve  these  problems  diverse 
numerical  methods  have  been  reported  which  include 
finite  difference  [4],  finite  element  [5],  and  boundary 
integral  methods  (BIM)  [6-8].  Using  these 


conventional  methods,  it  has  been  found  that  fine  grids 
and  a  large  number  of  elements  must  be  employed  to 
get  satisfactory  accuracy  [3].  This  requires  large 
computer  core  storage,  and  more  computational  time 
especially  for  the  iteration  scheme  of  the  nonlinear 
Poisson’s  problem  where  the  value  at  each  grid  point 
needs  to  be  updated  at  each  step  of  the  iteration. 
Further,  the  BIM  formulations  are  in  most  cases 
limited  to  homogeneous  Helmholtz  equation  and  tied 
closely  to  the  particular  problem  at  hand  [6].  In  this 
paper,  a  simple  approach  to  solve  the  homogeneous 
and  nonhomogeneous  Helmholtz  equations  is 
proposed.  The  technique  is  based  on  the  computation 
of  Laplacian  potential  by  the  method  of  moments 
(MoM)  [9],  without  resorting  to  different  formulations 
using  Hankel  functions  as  it  is  commonly  done  in  BIM 
[10].  Besides  its  generality  to  solve  Laplace’s, 
Poisson’s,  and  Helmholtz’s  equations  in  one  single 
code  implementation,  the  present  method  will 
considerably  reduce  the  number  of  domain  grids 
compared  to  the  finite  difference  methods  and  does 
not  require  any  interpolation.  The  accuracy  of  the 
MOM  solution  from  this  new  formulation  will  be 
compared  with  the  solutions  of  fine  difference  method 
(FDM)  and  finite  element  method  (FEM),  using  the 
ELLPACK  implementation  [4]. 

II.  MATHEMATICAL  FORMULATION 

Consider  the  following  two-dimensional  elliptic 
equation  for  a  smooth  function  defined  in  a  2-D 
region  defined  by  91  which  is  bounded  by  a  contour  C 
so  that 

Vlx¥(x,y)+'k(x,y)x¥(x,y)=F(x,y )  (1) 
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where  X  and  F  are  known  functions  in  the  domain  9L 
The  general  form  of  (1)  includes,  as  specializations, 
the  following  cases: 

1)  Laplace’s  equation,  with  X  =  0  and  F  =  0 

2)  Poisson’s  equation,  with  X  =  0  and  F *  0. 

3)  Helmholtz’s  equation,  with  X  *  0  and  F  *  0. 
So  only  one  general  formulation  addresses  the  solution 
to  all  of  the  three  cases  above.  In  addition  we  illustrate 
how  to  use  the  frequency  independent  solution  of  the 
Laplace’s  equation  to  solve  the  general  Helmholtz 
equation.  On  the  contour  C  the  boundary  condition  can 
be  of  Dirichlet,  Neumann,  or  mixed  type,  as  given  by 
the  general  form 


(x'F+p— =y  (2) 

on 

where  a,  (3,  y,  are  known  spatial  functions  and 
represents  the  normal  derivative.  It  is  implied  that  a 
consistent  set  of  boundary  conditions  has  been  chosen 
for  the  problem. 

The  proposed  scheme  to  solve  the  given  boundary 
value  problem  starts  by  assuming  the  term  X*¥  to  be  a 
known  function,  and  including  it  with  the  given 
excitation  function  F,  and  thereby  reducing  (1)  to  the 
familiar  Poisson’s  equation 

v2my)=-©(^y)  (3) 

with 

0(*,>O  =  A(.*,;y)¥(x,;y)-F(.*,;y).  (4) 


0M> ')  = 


—  JJ©(/,  y')(n 
2-71  „ 


K 


yj(x-x)2  +(y-/)2 


px'dy' 


(8) 


where  the  spatial  sets  ( x ,  y )  and  (x\  y')  denote  the 
spatial  coordinates  of  the  field  and  source  coordinates, 
respectively,  and  K  is  an  arbitrary  constant.  Here, 
£n(  K  )  represents  the  value  of  the  scalar  potential 
at  infinity  for  the  two  dimensional  case.  For  the  3D 
case,  this  term  does  not  exist,  as  the  potential  at 
infinity  is  zero.  Hence  the  Green’s  function  is  simply 
1/R  instead  of  the  natural  log  function.  The  value  of 
the  parameter  K  is  chosen  to  be  100  for  the  2-D 
problem  as  the  reference  potential  at  infinity  is  not 
zero  but  finite. 

An  expression  similar  to  (8)  can  be  derived  to 
approximate  the  Laplacian  potential  <)>/,.  The  potential 
<})*  can  be  assumed  to  be  produced  by  some  equivalent 
sources  consisting  of  electrical  charges,  o  located  on 
the  contour  C  [12].  Then  the  potential  ^  at  any  point 
(x,  y)  can  be  obtained  from  oQt,  y)  using  the  following 
integral  [9] 


i 


in 


1  a{x,y’)t n 


K 


d(x-x)2+(y-y)2  J 


dl' 


(9) 


The  solution  to  Poisson’s  equation  in  (3)  can  be 
expressed  as 


^=(/>h+t  (5) 

where  is  the  solution  to  the  homogeneous  Poisson’s 
equation  (Laplace’s  equation) 


where  /'  is  the  arc  length  on  the  contour  C. 

The  boundary  condition  of  the  homogeneous 
potential  <|)/f  is  obtained  from  (2)  and  (5)  as 


a0h+  P 


dn 


(  dd)  ^ 


(10) 


o 

11 

> 

(6) 

and  is  the  particular  integral,  i.e., 

VVP  =  -8(*,  y) 

(7) 

Here  we  use  the  particular  solution  of  the  Poison’s 
equation  given  by  [12] 


It  can  be  seen  that  (6)  along  with  the  boundary 
condition  of  (10),  constitute  a  similar  boundary- value 
problem.  A  similar  problem  was  addressed  in  [9] 
almost  thirty  years  ago. 

This  completes  the  formulation  of  the  problem. 
The  characteristic  features  of  this  formulation  are: 

a.  The  frequency  term  appears  in  an  explicit 
form. 

b.  The  Green’s  function  and  the  unknowns  are 
independent  of  frequency.  The  Green’s 
function  is  in  fact  the  static  Green’s 
function. 
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Next  we  illustrate  how  to  solve  these  coupled 
equations  both  for  the  homogeneous  part  and  the 
particular  integral. 

III.  NUMERICAL  SOLUTION  USING 
THE  METHOD  OF  MOMENTS  (MOM) 

In  order  to  solve  the  above  equations  given 
by  (8)-(10)  we  employ  the  Method  of  Moments. 

To  evaluate  the  above  integral  in  (8)  we 
divide  the  domain  9t  into  N  sub-regions.  The 
midpoint  coordinates  of  each  of  the  sub-regions  A% 
are  denoted  by  (x2i,  y2/)-  Then  the  potential  §p  at  the 
field  point  (x,  y )  is  approximated  by 


c,(x,y)  = 


Using  (4),  (11)  and  (14),  the  Helmholtz  potential  at 
an  arbitrary  field  point  (x,  y)  can  now  be  expressed  as 
'¥(x,y)  = 

M  N 

<  =  1  i  =  1 

(16) 


<t>P  (x,  y)=^  ©(*2, .  y2i  M  (*.  y)  (ii) 


where 

at{x,y)- 


r  JJ<nf  i — tt — fr 

lK  1%  J(x-x)  +(y-y) 


dx'dy 


Hence  it  is  assumed  that  0(;t,  y)  is  constant  within 
each  sub-region  A%  and  is  equal  to  the  value  0 

fej2/). 

Next  we  show  how  to  numerically  evaluate  (9). 
Here,  we  take  recourse  to  MOM.  We  expand  the 
surface  unknown  a  using  a  pulse  expansion.  For 
purposes  of  illustration  and  simplicity  let  us  choose 
point  matching.  Through  the  use  of  a  pulse-expansion 
and  point-matching  techniques  we  will  solve  the 
present  problem  given  by  (9).  Furthermore,  if  the 
contour  C  is  segmented  by  M  straight  lines  of  length 
AC,  between  points  i  and  i  +  1,  then  a  can  be 
represented  by  the  step  approximation 

M 

*=Z<rf,(0  (13> 

i  =  1 


where  for  simplicity  we  have  used  the  following 
abbreviations  in  (4): 

h  =  k  (*2 h  yid 
F,  =  F  (x2h  yii) 

%=  ^(X2iy2i) 

The  unknown  function  ¥(x,  y)  in  (16)  of  Helmholtz 
equation  can  now  be  evaluated  once  the  unknown 
terms,  g2  and  xFf,  are  determined.  Next  a  system  of 
two  matrix  equations  is  derived  and  solved  for  the 
unknowns  a,  and  The  first  matrix  equation  of  this 
system  is  readily  obtained  by  satisfying  (16)  at  the 
midpoints  (x2ji  y2j)  of  each  of  the  N  sub-region  A9ti. 
Using  matrix  notation,  we  obtain 

W®[Pjt][°il+[^][W-^]  (17) 

where  p;f  =  Ci(x2j,  y2j),  and  q}i  =  afyty  y2j). 

The  second  matrix  equation  is  now  obtained  by 
enforcing  the  boundary  condition  of  (10)  as  follows. 

We  define  fj  -  (Xjt  yj),  j  =  1,  2,  M,  to  be  the 

midpoints  of  AC;-.  The  boundary  conditions  are 
enforced  at  each  fj .  Substitution  of  (1 1)  and  (14)  into 

(10)  gives  the  following  set  of  equations 


where  P ,(/)  is  a  pulse  function  equal  to  1  on  A C,  and 
zero  elsewhere  and  o,  is  its  unknown  amplitude. 
Substituting  (13)  into  (9),  we  obtain  an  approximation 
for  4 )h 

M 

$k{x*y)  =  ^dOffi{xty)  (14) 

i  =  1 

where 


(18) 


(19) 
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h  = 


aai  +  p 


dal 

dn 


(20) 


Ax,y)=fj 


Equation  (18)  can  now  be  written  in  matrix  notations 
as 


<2|> 


illustrated  through  numerical  examples  in  the 
next  section. 

5)  Once  is  determined  at  the  grid  points, 
using  Gaussian  elimination  for  instance  the 
potential  at  any  other  point  is  obtained  using 
ordinary  matrix  multiplications  as  shown  by 
(17).  No  interpolations  are  needed. 

Next,  we  illustrate  these  points  through  some 
numerical  examples. 


Observe  that  (17)  and  (21)  form  a  system  of  two 
equations  in  two  unknowns  o,  and  4*,,  from  which 
they  can  be  solved.  We  use  (21)  to  obtain  an 
expression  for  a  and  then  substitute  it  into  (17),  and 
after  simple  matrix  manipulations  we  obtain  the 
following  equation  for  x¥l 

=  W  (22) 

where 

(23) 

and  [/]  denotes  the  N  x  N  identity  matrix. 

From  the  last  three  equations  some  general 
comments  can  be  made: 

1)  Analytic  expressions  can  easily  be  derived 
for  the  evaluation  of  all  the  terms  of  the 
matrices  [A]  and  [/?]  that  include  the  integrals 
(12)  and  (15).  This  is  true  at  least  for  the  2-D 
case. 

2)  The  matrix  elements  [vv;/],  which  are  a  part  of 
the  system  matrix  of  the  moment  method,  are 
similar  to  the  matrix  elements  obtained  in  [9]. 
For  different  Helmholtz  problems  with 
different  boundary  conditions,  [vty]  remains 
unchanged.  This  observation  is  important  in 
an  iteration  scheme  where  the  boundary 
conditions  are  kept  the  same. 

3)  Problems  with  multiple  right-hand  sides  are 
solvable  with  minimum  additional 
computation  time  since  the  function  F  in  (17) 
appears  only  as  a  term  of  matrix  [/?]  in  (24). 

4)  The  domain  and  the  contour  discretization 
schemes  are  totally  independent.  This  feature 
can  offer  the  flexibility  to  handle  boundary 
discontinuities  without  the  need  for  excessive 
domain  grid  generations.  This  will  be 


III.  STATIC  APPLICATIONS 

Example  1 — Water  Wave  Propagation :  As  a  first 
example,  we  consider  the  problem  of  water  wave 
propagation  in  a  rectangular  basin  100m  x  100m  [3]. 
Denoting  by  W  the  water  evaluation,  then  the  wave 
propagation  is  governed  by  (1)  with  F  =  0,  and  X  =  k2 
where  k  is  the  wave  number.  The  boundary  conditions 
used  are  displayed  in  Figure  1.  The  solution  by  MOM 
is  compared  to  the  solutions  obtained  by  finite 
difference  method  (FDM)  and  finite  element  method 
(FEM)  for  various  mesh  sizes  inside  the  domain.  The 
minimum  number  of  nodes  required  for  each  method 
to  converge  to  the  exact  solution  at  any  arbitrary  point 
is  optimized.  Figure  2  compares  the  solutions 
computed  along  the  line  x  =  90  by  the  Method  of 
Moments  utilizing  (lxl),  (2x2)  and  (4x4)  grids. 
Figure  3  shows  the  solutions  obtained  by  the  finite 
element  method  utilizing  (5x5),  (10x10),  and  (20x20) 
grids.  Figure  4  shows  the  results  obtained  by  the  finite 
difference  method  utilizing  (5x5),  (15x15),  and 
(35x35)  grids.  In  all  these  figures,  the  solid  line 
represents  the  exact  solution.  Referring  to  these 
figures,  one  can  easily  observe  that  a  considerably 
smaller  number  of  nodes  is  employed  with  the  present 
method  than  with  the  conventional  methods.  The 
minimum  number  of  the  domain  nodes  required  to 
obtain  a  satisfactory  accuracy  using  MOM  for  the 
present  method  is  16  (4x4),  compared  to  400  (20x20) 
for  FEM  or  up  to  1000  nodes  for  FDM. 

Example  2 — MOST  Modeling :  The  nonlinear 
Poisson’s  equation  plays  a  key  role  in  numerical 
modeling  of  semi-conductor  devices.  Many  important 
characteristics  of  VLSI  devices  can  be  extracted  from 
the  solution  of  Poisson’s  equation.  The  most  common 
approach  to  the  numerical  solution  of  the  nonlinear 
Poisson’s  equation  is  based  on  the  application  of 
Newton’s  method  to  simultaneous  discretized 
equations  [15].  This  approach  often  requires  large 
storage  especially  for  fine  meshes,  as  is  the  case  for 
two-dimensional  modeling  of  the  MOST  [14]. 
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*=ft(x) 
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where  f  (x)  = 


$=fb(x) 

cos  (ax)  if  20  <  x  <  80 
-2  if  x  <  80 


1  if  x>80 

fs(x)  =  sin  (ay),  fb(x)  =  cos(ax),  and  a=2rc 

Figure  1.  Physical  structure  and  the  associated  boundary  conditions 
associated  with  Example  1. 


Figure  2.  Solution  obtained  by  the  present  method. 


Figure  3.  The  FEM  solution  for  Example  1. 


In  this  section  the  MOM  is  applied  to  solve  the 
nonlinear  Poisson’s  equation  that  arises  in  the  MOST 
modeling.  We  consider  the  MOST  structure  of  Figure 
5  made  on  a  p-type  substrate  with  doping  NA.  Under 
the  low  current  approximation  the  potential  O  is 
governed  by  the  Poisson’s  equation  [16] 


a2y  82y  _  nf 

dx2  By2  l}D 


i  N 

e  —e  — 


A 

ns 


(25) 


where  ¥  =  $>IVT  is  the  normalized  potential,  VT  is  the 
thermal  potential,  nt  is  the  intrinsic  carrier 
concentration,  and  LD  is  the  Debye  length. 


<b  =  V  -  V 

^  VG  Vfb 


j 

k  t 
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ox 

ox 

_ 3 
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<f>  ss  V  -  V 

^  VD  v  si 


II 

<8 


$  =o 

Figure  5.  The  physical  structure  and  the  associated  boundary 
condition  for  MOST  modeling. 

The  boundary  conditions  adopted  along  the  edges 
of  the  device  are  the  same  as  those  used  by  [16],  [12], 
and  are  displayed  in  Figure  5.  On  the  oxide- 
semiconductor  interface  the  following  boundary 
condition  is  assumed 


f  r  9* 

hx  dy 


(26) 
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where  Eox ,  and  Es  are,  respectively,  the  oxide  and 
silicon  permittivities,  VFB  is  the  flat  band  voltage,  and 
tox  is  the  oxide  thickness. 

To  solve  the  above  elliptic  problem  we  proceed 
first  by  dividing  the  boundary  into  M  segments.  Finer 
segments  are  used  on  the  top  edge  of  the  device  to 
handle  its  boundary  discontinuities.  Using  the 
boundary  conditions,  and  the  technique  described  in 
[9],  matrix  [vty-]  is  computed,  then  inverted  and  stored. 
Independently,  the  base  region  is  also  divided  into  N 
finite  elements  or  cells.  We  seek  to  determine  the 
electric  potential  at  each  midpoint  of  the  cells  using 
(25)  which  is  nonlinear.  To  solve  it  we  set  up  an 
iterative  procedure,  based  on  Newton’s  linearization 
[15].  At  the  £-th  iteration  we  replace  the  right-hand 
side  of  (25)  by  its  Fourier  expansion  about  ¥*.  Then 
the  following  Helmholtz’s  equation  is  obtained 

02\j/(*+l)  02vp(*+l) 

dx2  dy2 

F(W(k)) 

(27) 

where  F  and  G  are  given  by 

F(4/W)  = 

^r((¥k)  +\)e-'¥a)  +(¥k)  -l)^ 

ld  V.v  ni  J 

(28) 


(*+D 


+  e 


«(*)  1 


(29) 


The  iteration  scheme  starts  by  taking  some  initial 
guess  value  for  ^(0)  so  that  (27)  can  be  solved  for  the 
first  approximation  *P(1).  Then  4/(l)  is  used  to  find  the 
second  approximation  *P(2).  The  procedure  is  repeated 

until  the  norm  II  —  xft(k)  |j  -s  jess  tjian  a 


desired  tolerance.  However,  this  iteration  scheme 
often  diverges  [13],  [15],  and  some  damping  factor 
was  found  necessary  to  improve  the  convergence  of 
the  iterative  solution.  At  the  beginning  of  the  k- th 
iteration  step,  the  following  formula  was  used  for  all 
inner  nodes 


The  accuracy  of  the  solution  obtained  by  the 
present  technique  is  demonstrated  by  comparing  it 
with  the  FDM  (5-point)  solution.  For  both  methods, 
MOM  and  FDM,  we  let  the  values  of  the  electric 
potential  to  be  updated  at  each  mesh  point  by  means  of 
an  explicit  formula,  that  is,  without  the  solution  of 
simultaneous  algebraic  equations  [15].  We  assign  = 
0,  as  an  initial  guess  for  all  the  inner  nodes.  Once  the 
convergence  is  attained  for  the  domain  nodes,  then  the 
electric  potential  at  any  other  point  in  the  device  is 
determined  by  a  matrix  multiplication  as  shown  in 
(17),  and  without  the  need  of  any  interpolation.  For 
numerical  computations,  the  following  data  were  used: 
The  thickness  of  the  oxide  layer  was  tox  =  0.5  (im,  the 
flat  band  voltage  is  VFB  =  -1  V  and  the  doping  profile 
is  assumed  to  be  uniform  with  NA  =  1018  cm"3,  = 

1.5xl0,n  cm-3,  thermal  potential,  Vr  =  0.0258  V,  and 
R  =  0.1.  The  results  of  the  computations  are  shown  in 
Figure  6,  where  the  distribution  of  the  electric 
potential  at  the  thermal  equilibrium  is  plotted  along 
different  lines  parallel  to  the  Jt-axis.  The  solid  lines 
represent  the  solution  obtained  using  the  FDM  and  (o) 
symbol  is  reserved  for  the  solutions  obtained  by  the 
present  method.  Close  agreements  can  be  observed 
between  the  two  methods.  However,  for  FDM  721 
nonuniform  mesh  points  are  employed  to  reduce  the 
total  number  of  nodes  {for  uniform  meshes  over  4900 
(70x70)  nodes  ought  to  be  used}.  Finer  meshes  were 
chosen  in  the  depletion  region  and  near  the  junctions 
to  reproduce  accurately  the  fast  variation  of  the 
electric  potential  [16]  and  the  surface  discontinuities 
of  the  potential;  whereas  for  MOM,  the  number  of  the 
uniform  meshes  was  529  (23x23)  or  400  for 
nonuniform  cells.  The  number  of  iterations  required 
for  the  convergence  was  5 1  with  the  present  technique 
as  compared  to  108  iterations  with  FDM  to  reach  the 
point  at  which  the  absolute  maximum  between  two 
subsequent  iterations  was  less  than  0.05  (tolerance). 

Next  we  apply  this  method  to  the  solution  of  the 
cutoff-frequencies  of  various  waveguides. 


*  *0  Ol  0.2.  n.i  (1,4  o.S  TT*  7C*  . Vi/i . j 


tlHano.  (orr») 


\p(*+i)  =  (j  _  z?)^^  “  o  -j. 
where  R(<  1)  is  the  relaxation  factor. 


(30) 


Figure  6:  The  distribution  of  the  electric  potentials  along  the  lines 
parallel  to  the  x-axis. 
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IV.  SOLUTION  OF  THE  GENERAL 
HELMHOLTZ  EQUATION  FOR 
HOMOGENEOUSLY  FILLED 
WAVEGUIDES 

In  waveguides,  solution  of  the  Helmholtz 
equation  determines  the  electromagnetic  field 
configuration  within  the  guides.  It  is  convenient  to 
divide  the  possible  field  configurations  within  the 
waveguides  into  two  sets,  namely  TM  waves  and  TE 
waves,  each  of  which  is  governed  by  similar 
Helmholtz  equations. 

If  we  consider  a  waveguide  in  which  the 
direction  of  propagation  of  the  wave  is  along  the  z- 
direction,  then  the  Helmholtz  equations  are  as  follows. 

TM,  Case  (Hz  =  0): 

V2£z  (x,  y )  +  (g)2/ie  -  K2z)Ez  (. x ,  y)  =  0  (31) 

with  the  appropriate  boundary  conditions  Ez  =  0,  on 
the  conductor  walls  and  for  the 

TE,  Case  ( Ez  s  0): 

V2Hz  (x,  y)  +  (a)2  fi£  -Kz)Hz(x,y)  =  0  (32) 

with  appropriate  boundary  conditions  dHJdn  =  0,on 
the  conductor  walls  where  dHJdn  represents  the 
normal  derivative.  Here, 

Ez  z-component  of  the  electric  field; 

Hz  z-component  of  the  magnetic  field; 

co  angular  frequency  =  2 nfi 

f  frequency  of  interest; 

(i  permeability  of  the  homogeneous 

medium; 

£  permittivity  of  the  homogeneous 

medium; 

Kz  propagation  constant  in  the  z- 

direction. 


Comparing  (31)  and  (32)  with  (1),  and  also 
the  boundary  conditions  of  the  TM,  and  TE,  cases  with 
those  of  the  general  equations,  we  can  draw  the 
following  analogies  as  outlined  in  Table  I. 

By  examining  (32)  that 

[«]=[pJ,][w,r’[j'j 

it  can  be  inferred  that  [B]  =  0,  since  F  =  0  and  y  =  0  for 
TMZ  and  TEZ  cases. 

In  the  case  of  TMZ,  (22)  reduces  to  the  form 

[A]  •  [Ezi]  =  0  (33) 

In  the  case  of  TEZ,  (22)  reduces  to  the  form 

[A]-[ffZI-]  =  0  (34) 

Here  again,  Ezi  and  Hzi  refer  to  the  values  of  Ez  and  Hz 
at  the  midpoints  of  the  subregions  of  the  discretized 
waveguide  cross  section. 

For  (33)  and  (34),  nontrivial  solutions  exist 
for  [£,/]  and  [ Hzi ]  only  if  the  matrix  [A]  is  singular. 
The  condition  for  nontrivial  (i.e.,  nonzero)  solutions  to 
exist  for  [Ezi]  and  [Hzi]  it  is  essential  that 

det[A]  =  0  (35) 

where  det[A]  stands  for  determinant  of  [A],  We  know 
from  (23)  that 

[A]  =  ([p;i][w/J-1[&.i]-[9..])[^]+[ /] 

and  we  also  know  that  for  the  cases  of  the  TM,  and 
TE,  waves 

A  =  G)2 /HE  —  K2 .  (36) 


Table  I 

Special  cases  of  the  General  Helmholtz  Equation 


General  Equation  (3) 

TM,  Equation  (1) 

TE,  Equation  (2) 

V2'F+X4/=F 

V2Fz  +( £02|X  e  —  Kz  )-Ez=  0 

V2HZ+(0)2|4  8  -K2z)Hz=  0 

Hz 

X 

to2ji  8  -  K2 

C02|X  e  -K\ 

F 

0 

0 

t*1 

M 

II 

o 

dHz/dn=0 

7 

0 

0 
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Hence,  given  the  frequency  at  which  the  Helmholtz 
equation  is  to  be  solved,  det  [A]  would  be  a  function  of 
Kz,  the  roots  of  which  will  provide  the  values  of  Kz  for 
which  det  [A]  =  0.  Once  these  Kz  values  are  known, 
the  eigenvector  of  [A]  corresponding  to  the  minimum 
eigenvalue  gives  the  nontrivial  solutions  for  [Zy  and 
[Hzi]  in  case  of  TM~  and  TE.  cases,  respectively. 

Once  [Ezi\  and  [Hzi]  are  determined  at  the  grid 
points,  using  Gaussian  elimination  for  instance,  the 
values  of  Ez  and  Hz  at  any  other  point  can  be  obtained 
using  ordinary  matrix  multiplications,  as  explained  in 

tl]. 


IVa.  Calculation  of  Propagation  Constants  for 
Different  Waveguide  Modes 

It  is  evident  that  for  the  existence  of  nontrivial 
(nonzero)  solutions  for  [Ezi]  and  [Hzi]>  it  is  necessary 
that  (35)  be  satisfied.  Let  us  define  a  matrix  [Z]  such 
that 


[Z]  =  ([Pji][w, (37) 
Hence,  (35)  becomes 


Therefore,  the  propagation  constants  of  different 
modes  in  the  waveguide  are  given  by  the  following: 

For  (ATi)2 >0, 


*:=  CoV  8  + 


For  (K{ )2<0, 


EV{Z] 


propagating  modes 

(42) 


K'.=j  C02|i£  |2|  nonpropagating  modes 

'  EV; 


(43) 


Results  of  propagation  constants  of  various  modes  in  a 
rectangular  waveguide  computed  by  this  method  are 
shown  next. 


IV  b.  Calculation  of  Cutoff  Frequencies  for 
Different  Waveguide  Modes 

The  cutoff  frequencies  for  the  various  propagating 
modes  in  the  waveguide  are  given  by 


det  ([Z]  [A,,-]  +  [/])  =  0 

which  can  be  rewritten  as 


( 

f-0 

\ 

det 

[Z]- 

[/] 

V 

l 

J 

(38) 


(39) 


Equation  (39)  is  similar  to  the  characteristic  equation 
of  matrix  [Z],  with  its  eigenvalues  given  by  -lA*. 

Knowing  that  above  =C02|I  £  -  K\  for  TM,  and 
TE,  cases,  it  can  be  concluded  that 

ll  = _ ! _ -  EV. |Z> 

4  {K[)2-(02fl£  '  ’  (40) 

1  =  1, 2,...,N. 


where  Klz  is  the  propagation  constant  of  the  ith  mode 
and  EVj-Z^  =  /th  eigenvalue  of  [Z]. 


Equation  (40)  can  be  rearranged  as 


(Ki)2  =0)z/i£  + 


1 


EV}Z]  ' 


(41) 


f!  =  V(co2|i£  -  (K[)2)  (44) 


where  flc  =  cutoff  frequency  of  the  /th  mode.  Here, 
V  =  velocity  of  light  in  the  homogeneous  medium 

=  1  /  ^/ji  £  and  it  can  be  deduced  from  (41)  that 


co2  fje-i^Kiy 
i  =  1,2,.. 


-1 

“  EV}Z]  ’ 
.,N 


(45) 


Using  (45)  in  (44),  we  find  the  cutoff  frequencies  for 
the  first  N  propagating  modes  as 


f  = 


V 


-1 


2xyEV,lzr 


(46) 


The  cutoff  wave  number  klc  of  the  (th  mode  can  be 
calculated  from  the  cutoff  frequency  using  the  relation 

k‘c=^E,  1  =  1, 2,...,N 


V 


This  method  thereby  provides  a  straightforward 
approach  to  find  the  cutoff  frequencies  (and,  hence, 
cutoff  wavenumbers)  of  any  waveguide  structure 
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without  resorting  to  scanning  over  a  wide  range  of 
frequencies,  as  is  done  in  the  Ritz-Galerkin  and 
surface  integral-equation  methods. 

V  a.  Results:  Rectangular  Waveguide 

Consider  a  rectangular  waveguide.  For  the 
waveguide  in  Figure  7,  the  region  was  divided  into 
100  subregions  and  the  boundary  was  discretized  into 
96  subcontours.  The  maximum  matrix  size  involved 
in  the  computations  was  100  x  100.  Results  have  been 
displayed  in  Table  II  for  the  cutoff  wave  numbers  of 
the  first  eight  TM/TEZ  modes.  The  computational 
time  involved  in  finding  the  cutoff  wavenumbers  of 
the  first  100  modes  on  a  Sun  SPARC  10  workstation 
was  16  s. 


* 


t 


4  cm 

Figure  7.  A  rectangular  waveguide. 


Table  II 


Cutoff  wavenumbers  for  air-filled  rectangular  wave 

guide 

Mode 

No. 

Mode 

ke  actual 
(rad/cm) 

K 

computed 

(rad/cm) 

Diff. 

% 

1-0 

TEZ 

0.7857 

0.7921 

0.81 

0-1 

TE, 

1.0476 

1.0536 

0.57 

1-1 

TEZ, 

tm7 

1.3095 

1.3239 

1.00 

2-0 

te7 

1.5714 

1.5827 

0.72 

2-1 

TE,, 

TM. 

1.8886 

1.9108 

1.10 

0-2 

te7 

2.0952 

2.1095 

0.68 

1-2 

TMZ> 

TM, 

2.2377 

2.2610 

1.00 

3-0 

te7 

2.3571 

2.3896 

1.30 

Vb.  Single-Ridge  Waveguide 

A  single  ridge  waveguide  is  a  popular  means  of 
getting  higher  bandwidth.  The  first  four  TMz  and  TEz 
mode  cutoff  wavenumbers  were  computed  for  the 
single-ridge  hollow  waveguide  shown  in  Figure  8. 
Results  have  been  displayed  in  Table  III  and  compared 
with  published  data.  For  the  waveguide,  the  region 


was  divided  into  96  subregions  and  the  boundary  was 
discretized  into  112  subcountours.  The  maximum 
matrix  size  involved  in  computations  was  96  x  96. 
The  computational  time  involved  in  finding  the  cutoff 
wavenumbers  of  the  first  96  modes  in  each  case  on  a 
Sun  SPARC  10  workstations  was  18  s. 


Table  III 

Cutoff  wavenumbers  for  air-filled  single-ridge  waveguide 


Mode 

No. 

Mode 

kc 

published 

(rad/cm) 

kc 

computed 

(rad/cm) 

Diff. 

% 

i. 

TM, 

12.16405 

12.2338 

0.57 

2. 

tm7 

12.293817 

12,4106 

0.95 

3. 

™7 

13.996417 

14.2152 

1.56 

4. 

tm7 

15.587117 

15.8221 

1.50 

5. 

te7 

2.2566s 

2.2688 

0.54 

6. 

te7 

4.943617 

5.0149 

1.44 

7. 

te7 

6.518917 

6.6289 

1.68 

8. 

te7 

7.564217 

7.7097 

1.92 

VI.  CONCLUSION 

An  efficient  technique  based  on  MoM  formulation 
for  solving  a  general  Helmholtz  equation  starting  from 
Laplace’s  equation  is  presented.  The  main  feature  of 
this  new  formulation  is  that  the  boundary  conditions 
are  satisfied  independent  of  the  discretizations  of  the 
regions  and  the  nodes.  This  feature  was  found 
especially  useful  when  the  boundary  conditions  have 
discontinuities.  Considerable  reduction  in  the  domain 
grids  is  realized  with  the  present  method  as  compared 
to  the  conventional  methods  such  as  finite  difference 
method  or  the  finite  element  methods.  In  addition,  one 
need  not  use  a  frequency  dependent  Green’s  function 
which  can  reduce  the  computational  cost  significantly. 
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Abstract  -  In  this  paper  we  present  a  simple 
modification  to  the  traditional  Finite  Difference 
Time  Domain  (FDTD)  method  for  treating  material 
cells.  The  Yee  cell  is  split  into  8  smaller  material 
subcells  so  that  each  E  and  H  field  point  is 
considered  to  be  located  at  the  crosspoint  of  8  cells 
with  differing  material  properties.  Thus  there  is  no 
longer  an  overlap  of  the  material  cells  associated 
with  the  components  of  the  E  and  H  fields.  The  8 
materia]  properties  are  averaged  at  each  crosspoint. 
Since  the  averaging  is  done  outside  the  time 
marching  loop  there  is  little  increase  in  the  total 
computational  time.  Numerical  results  are  shown  for 
a  sinusoidal  plane  wave  scattering  from  a  dielectric 
sphere.  These  results  are  compared  with  the  exact 
Mie  solution  and  the  traditional  material  cell  method 
along  different  cuts  through  the  sphere.  The  splitting 
and  averaging  is  shown  to  give  improved  amplitude 
accuracy  in  the  vicinity  of  the  sphere.  Improvement 
is  also  observed  at  planar  interfaces  angled  with 
respect  to  the  grid.  An  additional  benefit  of  this 
subcell  formulation  is  that  objects  may  be  modeled 
with  twice  the  geometrical  resolution  without 
increasing  the  size  of  the  staggered  grid. 

Index  Terms-  FDTD,  Subcell,  Split  cell,  Effective 
Dielectric  Constant 


I.  INTRODUCTION 

The  finite-difference  time  domain  (FDTD) 
technique  is  a  widely  used  method  for  computing 
electromagnetic  scattering.  The  formulation  is  based  on 
the  ideas  put  forth  in  [1]  and  involve  computing  the 
electric  and  magnetic  (E  and  H)  fields  on  a  staggered 
grid.  There  are  some  shortcomings  with  the  Yee  method 
when  applied  to  interfaces  between  materials.  First,  the 
material  cells  associated  with  the  field  components 
(staggered  in  space)  are  centered  around  each  field  point 
and  thus  overlap  in  space.  Second,  the  discontinuity  in 
permittivity  at  interfaces  between  different  media  causes 
accuracy  problems.  Third,  is  the  error  caused  by 
approximating  a  curved  boundary  by  a  staircased  one. 

Various  methods  have  been  proposed  to  improve 
the  accuracy  of  the  basic  FDTD  method  at  interfaces, 
mainly  by  reducing  the  staircasing  errors  caused  by 
angled  or  curved  boundaries.  For  example,  one  method 
uses  a  volume  weighted  average  of  €  on  each  side  of  a 


material  interface  [2].  An  effective  dielectric  constant  is 
computed  as  follows: 

^eff  =  (^1^1  ^2£2  )/^cell  >  0) 

where  V}  is  the  volume  on  one  side  of  the  interface  and 

V2  is  the  volume  on  the  other  side  of  the  interface.  This 

method  requires  a  special  treatment  of  each  material 
interface.  Specifically  the  method  requires  identifying 
which  cells  are  intersected  by  the  boundary  and  where 
the  boundary  intersects  these  cells.  Dey  and  Mittra 
reported  that  the  results  using  this  approach  are 
significantly  more  accurate  than  using  a  staircase  to 
approximate  the  surface  of  a  dielectric  structure. 
Another  approach  uses  effective  dielectric  constants 
together  with  the  contour-path  integral  FDTD  (CFDTD) 
method  [3].  Recently  a  new  technique  was  developed 
that  utilizes  the  electric  field  components  along  the 
edges  of  the  cell  [4],  Rather  than  using  the  volume  of 
the  dielectric  occupying  a  cell  it  instead  uses  the  length 
of  the  dielectric  intersecting  the  cell  edges. 

Still  another  method  uses  dielectric  subgrid 
resolution  (DSR)  [5].  This  can  be  described  as 
partitioning  a  standard  cell  into  8  subcells,  each 
homogeneously  filled.  The  8  cells  can  have  different 
material  properties.  Generalized  constitutive  parameters 
are  obtained  by  treating  each  cell  as  an  equivalent 
lumped  circuit.  It  is  shown  by  numerical  example  that 
on  coarser  grids  the  simpler  averaging  method  gives 
more  accuracy  than  the  lumped  circuit  approach  across 
interfaces.  The  DSR  method  also  requires  greater 
memory  since  there  are  effective  dielectric  constants  for 
3  directions. 

The  practical  advantage  to  the  approach  presented 
in  this  paper  is  that  no  special  preprocessing  operations 
are  required  to  compute  the  intersections  of  cells  with 
interface  boundaries.  The  comparisons  in  this  paper  are 
on  the  basis  of  field  amplitudes  as  opposed  to  resonance 
frequencies  used  in  the  cited  articles. 

II.  SIMPLE  SUBCELL  MODEL 

The  method  presented  here  uses  the  standard 
staggered  grid  scheme  but  now  sub-divides  or  splits  the 
existing  cell  into  8  smaller  cells  as  shown  in  Fig.  1. 
Each  field  point  may  now  be  viewed  as  a  crosspoint 
surrounded  by  8  subcells  with  differing  parameters  as 
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opposed  to  being  at  the  center  of  a  material  cell.  The 
material  cells  no  longer  overlap  in  space  for  the 
staggered  field  points.  A  sphere,  for  example,  cutting 
through  the  standard  Yee  cell  still  is  approximated  with 
a  staircase  but  with  twice  the  number  of  cells.  Unlike 
some  previous  approaches  there  are  no  special 
calculations  related  to  the  specific  geometry  modeled. 

Applying  Ampere’s  Law  around  a  contour  Cx  one 
obtains: 


direction  as  shown  in  Fig.  2.  The  dielectric  sphere  has  a 
radius  of  4.5  cm.  In  this  work  we  consider  scattering  at 
two  frequencies  and  two  dielectric  constants.  Case  I  is  at 

a  frequency  of  2.5  GHz  with  £  —  4 £0  in  the  sphere  and 

Case  II  is  at  1.768  GHz  with£  =  $£0  in  the  sphere.  In 

each  case  the  wavelength  in  the  sphere  is  the  same.  The 
Mie  series  solution  for  the  fields  inside  and  outside  the 
sphere  is  compared  to  numerical  FDTD  results. 


(2) 

™  S,  CY 

From  Fig.  1,  assuming  that  Dz  is  constant  over  a 

surface  patch  S]  (center  x-y  plane  in  Fig.  1)  and  is 
equal  to  the  value  at  the  center  of  the  patch,  the  surface 
integral  becomes: 


J  D:dS  =  e(i,  j,  k)E:  |  ijk  Ax  Ay .  (3) 

Assume  also  that  the  value  of  the  H  field  along  each 
contour  leg  is  constant  and  equal  to  the  value  at  the 
midpoint  of  the  leg.  Since  the  material  cell  is  now  split, 
further  assume  that  at  the  cell  center  or 

crosspoint  of  the  8  material  cells,  that  the  effective 
dielectric  constant  is  the  average  of  the  8  surrounding 
material  cells  or : 


8 

e(i,j,k)  =  e{i,j,k)  =  i/8^£„ .  (4) 

n=\ 


Using  central-difference  approximations  for  the  time 
derivatives  (interleaved  in  time)  the  following  step- 
ahead  equation  results  -  valid  for  non-cubical  cells. 


J-V7+1 

ijjk 


KjJk 


A? 

+  ~ - 

£(i,j,k) 


ij-yw 


1  rr«+l/2 1 

Ay  *  I'W 


(5) 


Analogous  equations  can  be  developed  for  Ev  and  Er . 

y  X 


III.  BENCHMARK  PROBLEM  -  DIELECTRIC 
SPHERE 

The  problem  considered  is  a  plane  wave  scattering  from 
a  dielectric  sphere.  The  incident  wave  has  the  electric 
field  polarized  in  the  z  direction  and  travels  in  the  +y 


A.  Computational  Modeling 

Numerical  results  are  generated  using  a  3D  FDTD 
computational  code  that  implements  the  standard  Yee 
method  with  and  without  the  subcell  averaging.  The 
computational  grid  is  512  x  256  x  512.  At  the 
extremities  of  the  grid  Liao  second  order  absorbing 
boundary  conditions  (ABCS)  are  employed.  A  sine 
wave  source  is  fed  in  at  y=5.  The  incident  wave  is 
planar  over  1/2  of  the  grid’s  width  or  from  x=  128  to 
x=384,  with  a  Gaussian  taper  toward  the  absorbing 
boundaries  This  finite-width  waveform  approximates 
plane  wave  conditions  in  the  vicinity  of  the  sphere.  The 


spatial  discretization  Ax  is  set  to 


the 


dielectric  media)  for  both  Cases  I  and  II.  The  sphere  has 
a  radius  of  1 5  Ax  ;however,  with  the  split  cell  model, 
the  sphere  has  a  radius  of  30  subcells.  The  FDTD 
simulation  is  run  in  continuous  wave  (cw)  mode,  long 
enough  to  achieve  steady-state  conditions  over  the 
extent  of  the  grid.  The  amplitudes  of  the  vector 
components  of  the  E  field  are  saved  along  chosen  line 
cuts  (along  y)  though  the  sphere.  The  program,  which 
updates  large  3D  grids,  is  written  with  OpenMP 
directives  to  run  efficiently  on  an  SGI  Origin  parallel 
computer. 


B.  Analytical  Solution 

The  Mie  solution  for  a  plane  harmonic  wave 
scattering  from  a  dielectric  sphere  is  well  known  and 
may  be  found  in  several  references  [6,7].  This  exact 
solution  is  a  series  for  the  scattered  field  exterior  to  the 
sphere  and  the  transmitted  or  total  field  interior  to  the 
sphere.  The  incident  field  is  added  to  the  scattered  field 
to  obtain  the  total  field  outside  the  sphere.  The  resulting 
vector  components,  in  spherical  polar  coordinates,  are 
transformed  into  Cartesian  vector  components  along 
various  cuts  though  the  sphere.  This  permits  a  direct 
comparison  with  the  FDTD  solution.  Care  must  be 
taken  to  insure  that  the  FDTD  cuts  correspond  exactly 
in  space  with  the  analytical  ones. 


C.  Results  and  Comparison 

Fig.  3  shows  comparisons  for  Case  I  between  the 
Mie  solution  (exact),  the  8  subcell  or  split  cell  model, 
and  standard  cell  model  for  the  magnitude  of  Ez  along  a 
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straight  line  through  the  sphere.  For  cuts  parallel  to  the 
y-axis,  Ez(0,y,0)  in  Fig.  3(a)  and  Ez(8,y,0)  in  Fig.  3(b), 
the  standard  cell,  the  split  cell  and  exact  results  are  close 
and  difficult  to  distinguish.  However,  the  standard  cell 
model  exhibits  a  local  amplitude  overshoot  at  the  back 
of  the  sphere.  The  subcell  model  produces  more 
accurate  field  values  at  the  back  of  the  sphere. 

Fig.  4  shows  the  comparison  between  the  Mie 
solution,  the  subcell  model,  and  the  standard  model  for 
the  magnitude  of  Ey  Along  the  cut  Ey(0,y,l)  shown  in 
Fig.4(a),  the  standard  cell  model  gives  a  large 
undershoot  and  spatial  phase  error  near  the  back  of  the 
sphere,  where  there  is  a  discontinuity  in  Ey.  The  subcell 
model  does  not  exhibit  this  large  undershoot  or  spatial 
error.  The  standard  cell  model  also  overshoots  the  exact 
solution  near  the  front  of  the  sphere.  Fig.  4(b)  shows  the 
cut  Ey(0,y,8).  The  standard  cell  and  split  cell  are  closer 
in  amplitude  away  from  the  center  (z=0)  plane. 

Figs.  5  and  6  show  the  analogous  comparisons  for 
Case  II.  A  rather  large  overshoot  is  apparent  in  Fig.  5(b) 
near  the  back  of  the  sphere  in  the  z=8  plane  for  the 
standard  cell  case.  Fig.  6(a)  best  illustrates  the  standard 
cell  errors  and  shows  large  undershoots  at  both  the  front 
and  back  of  the  sphere,  where  there  are  jump 
discontinuities  in  Ey  An  overshoot  at  the  back  of  the 
sphere  for  the  standard  cell  method  is  also  apparent  in 
Fig.  6(b).  The  split  cell  model  agrees  better  with  the 
exact  at  both  the  front  and  back  of  the  sphere. 

To  compare  the  split  cell  and  averaging  method 
with  the  similar  dielectric  subgrid  resolution  (DSR) 
method  [5],  Case  II  was  computed  using  a  grid  spacing 
of  1/10  of  a  wavelength  in  the  sphere.  This  coarse  grid 
accentuates  any  differences  between  the  methods.  The 
Ey  component  across  the  interface  best  illustrates  these 
differences.  Fig.  7(a)  shows  that  the  DSR  method  tends 
to  underestimate  the  jump  discontinuities  as  compared 
to  the  simpler  averaging  method.  However,  at  1/20  of  a 
wavelength  in  the  sphere,  the  DSR  and  simple 
averaging  curves  are  both  converging  to  the  Mie  result, 
as  shown  in  Fig. 7(b). 

IV.  PLANAR  INTERFACES 

The  performance  of  the  split  cell  FDTD  method  is 
further  demonstrated  by  computing  the  reflection  and 
transmission  of  a  plane  wave  Gaussian  pulse  polarized 
in  the  z  direction  from  a  planar  interface.  The  grids  are 
chosen  to  be  128x64x128.  An  interface  tilted  at  14°  is 
simply  a  wedge  angled  from  (0,64,z)  to  (128,128,z). 
The  grid  spacing  is  based  on  the  2/3  down  power  point 
for  the  pulse.  At  this  frequency  the  grid  spacing  is 
chosen  to  be  1/10  of  a  wavelength  in  the  dielectric.  At 
smaller  grid  spacings  the  differences  between  the  split 
cell  and  standard  cell  become  minimal  for  a  single 
pulse.  For  1/10  of  a  wavelength  grid  spacing  in  the 
dielectric  a  second-order  accurate  FDTD  formulation  is 
used.  The  same  cases  are  repeated  with  a  coarser  grid 


spacing  of  1/7.5  of  a  wavelength  but  using  a  fourth- 
order  (second-order  in  time  and  fourth-order  in  space  or 
2,4  scheme)  FDTD  formulation. 

The  peak  amplitudes  of  the  transmitted  pulse  and 
reflected  pulse  are  compared  with  the  exact  amplitudes. 
In  order  to  capture  the  reflected  amplitude  the  incident 
pulse  is  subtracted  out  on  the  reflecting  side  of  the 
interface.  The  results  are  summarized  in  Tables  I-IV. 
The  percent  errors  are  shown  in  parentheses.  The  data 
show  that  the  split  cell  formulation  is  consistently  more 
accurate  for  interfaces  angled  at  14°.  The  cell  splitting 
and  averaging  is  particularly  useful  on  coarse  grids 
where  staircasing  errors  may  be  more  of  a  concern.  The 
split  cell  and  standard  cell  give  the  same  transmitted  and 
reflected  amplitudes  for  0°  (normal  incidence)  thus 
illustrating  that  the  averaging  only  has  an  effect  for 
interfaces  at  angles. 

V.  CONCLUSIONS 

The  split  cell  and  averaging  method  is  shown  to 
give  improved  accuracy  across  angled  interfaces  and 
along  curved  interfaces,  such  as  a  sphere.  The  method  is 
especially  useful  where  coarser  grids  are  required.  The 
technique  is  very  simple,  only  requiring  that  the 
material  cells  have  twice  the  resolution  of  the  field  cells. 
Since  the  averaging  of  the  dielectric  constants  is  done 
prior  to  the  time  stepping  of  the  fields  there  is  little 
added  computational  cost. 
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Fig.  1  Diagram  of  Yee  cell  with  8  subcells  along  with  the  position  of  the  E  and  H  vector  components. 


z 


Fig.  2  Diagram  of  incident  wave,  dielectric  sphere,  and  coordinate  system. 
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Figure  4a  and  4b  -  Magnitude  of  Ey  parallel  to  y-axis  through  (x=0,z=l)  and  (x=0,z=8).  Case  1,  f=2.5GHz  and  e=4. 
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Figures  5a  and  5b  -  Magnitude  of  Ez  parallel  to  y-axis  and  through  (x=0,z=0)  and  (x=0,  z=8).  Sphere  is  located 
between  y=l  13  and  y=143.  Case  II,  f=1.767  GHz  and  e=8. 
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Figures  6a  and  6b  -  Magnitude  of  Ey  parallel  to  y-axis  through  (x=0,z=0)  and  (x=0,z=8).  Sphere  is  located  between 
y=l  13  and  y=143.  Case  II,  f=l .767  GHz  and  e=8. 
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Figure  7a  -  Magnitude  of  Ey  parallel  to  y-axis  through  (x=0,z=l).  Compared  are  DSR  and  simple  averaging.  Case  II, 
f=  1.767  GHz,  and  e=8,  10  points/wavelength.  Figure  8b  is  same  as  8a  but  with  20  points/wavelength. 
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Table  I 

Amplitudes  of  Transmitted  and  Reflected  Waves 
Second-Order,  e=4.0,  14°,  A/10 


Split  Cell 

Standard  Cell 

Exact 

Transmitted 

Amplitude 

.660 

(.3%) 

.655 

(1.06%) 

.662 

Reflected 

Amplitude 

-.335 

(3.7%) 

-.344 

(6.5%) 

-.323 

Table  II 

Amplitudes  of  Transmitted  and  Reflected  Waves 
Second-Order,  £=8.0,  14°,  A./10 


Split  Cell 

Standard  Cell 

Exact 

Transmitted 

Amplitude 

.519 

(0.0%) 

.513 

(1.16%) 

.519 

Reflected 

Amplitude 

-All 

(1.9%) 

-.484 

(3.4%) 

-.468 

Table  III 

Amplitudes  of  Transmitted  and  Reflected  Waves 
Fourth-Order,  £=4.0, 14°,  A/7.5 


Split  Cell 

Standard  Cell 

Exact 

Transmitted 

Amplitude 

.660 

(0.3%) 

.643 

(2.8%) 

.662 

Reflected 

Amplitude 

-.325 

(.6%) 

-.351 

(8.7%) 

-.323 

Table  IV 

Amplitudes  of  Transmitted  and  Reflected  Waves 
Fourth-Order,  £=8.0,  14°,  A/7.5 


Split  Cell 

Standard  Cell 

Exact 

Transmitted 

Amplitude 

.518 

(0.3%) 

.502 

(3.3%) 

.519 

Reflected 

Amplitude 

-.468 

(0.0%) 

-.494 

(5.6%) 

-.468 
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Abstract 

The  aim  of  this  paper  is  to  outline  the  theory  for  calculating 
the  angular  and  radial  Mathieu  junctions  of  complex 
arguments.  These  functions  are  required  for  the  computation 
of  analytic  solutions  of  electromagnetic  scattering  by  lossy 
dielectric  elliptic  cylinders  and  waveguides.  The 
backscattered  echo  width  of  a  lossy  dielectric  elliptic 
cylinder  is  compared  with  the  special  case  of  lossy  circular 
and  weakly  lossy  elliptic  cylinders  and  excellent  agreement 
is  obtained  in  all  cases.  Tabulated  and  plotted  numerical 
results  of  typical  Mathieu  functions  are  presented. 

1.  Introduction 

Many  engineering  and  physics  applications  involve  the 
solution  of  structures  with  elliptical  geometries.  Analytic 
solutions  to  such  structures  require  the  computation  of 
angular  and  radial  Mathieu  functions  in  the  elliptical 
coordinate  system.  Examples  of  such  applications  are  the 
solution  of  elliptical  waveguide  problems  [1-5],  optical 
fibers  [6],  and  elliptical  loaded  horn  and  corrugated 
elliptical  horn  antennas  [2],  field  prediction  inside  lossy 
elliptic  cylinders  for  biological  body  modeling  [7],  and 
scattering  by  mutilayered  dielectric  elliptic  cylinders  [8-12]. 
Furthermore,  analytical  solutions  are  generally  needed  to 
check  the  accuracy  of  numerical  or  approximate  solutions. 
Over  the  years  many  algorithms  and  programs  have  been 
developed  to  address  the  problem  of  computing  Mathieu 
functions  of  integer  and  real  arguments  [4],  [13-22]. 
However,  none  of  these  publications  have  addressed  the 
more  general  problem  of  computing  the  angular  and  radial 
Mathieu  functions  of  complex  arguments.  Caorsi  et.  al. 
have  computed  the  solution  of  electromagnetic  scattering  by 
weakly  lossy  multilayer  elliptic  cylinders  using  first-order 
truncation  of  the  Taylor  expansion  of  each  Mathieu  function 
[23-24].  The  same  authors  have  presented  the  solution  of  the 
scattering  of  an  electromagnetic  wave  by  a  lossy  dielectric 
elliptic  cylinder  based  on  the  exact  Fourier  series  solution  in 


terms  of  the  Mathieu  functions  with  complex  argument  [24]. 
Recently,  Amendola  outlined  the  theory  for  the  computation 
of  Mathieu  functions  of  complex  order  by  applying  the 
Floquent  theory  [27].  Unfortunately,  there  are  no 
numerical  results  for  these  functions  given  in  the  literature. 
In  this  paper,  the  solution  presented  in  [15]  is  extended  to 
account  for  the  characteristic  values  (eigenvalues)  of 
Mathieu  functions  of  complex  argument.  Once  the 

characteristic  values  are  computed,  the  complex  Fourier 
coefficients  of  the  Mathieu  functions  can  be  obtained.  These 
coefficients  are  needed  to  compute  the  angular  and  radjal 
Mathieu  functions  of  complex  arguments  and  their 

derivatives  as  will  be  shown. 


2.  Theoretical  Background 


The  homogeneous  wave  equation  (Helmholtz)  in  elliptical 
coordinates  u,  v,  and  z  is  given  by 


_ 1 _ 

F2 (cosh2  w- cosh2  v) 


3V  dV 

+ 

[iU-l 

du2  8v2 

QJ 

N 

nj 

r  =  o  (1) 


The  elliptical  coordinate  system  ( u,v,z )  is  defined  in  terms  of 
the  Cartesian  coordinate  system  {x,y,z)  by 
x  =  Fcosh(w)cos(v)and  y  =  F sinh(w) sin(t') ,  where  F  is 
the  semi  focal  length  of  the  elliptical  cross  section. 
Assuming  no  variation  of  the  function  y/  along  the  z-axis. 


the  z-derivative  in  (1)  will  drop  out.  The  resulting 
differential  equation  can  be  separated  into  the  following 
differential  equations  in  terms  of  U(u)  and  V(v).  The 
circumferential  (ordinary)  Mathieu’s  differential  equation  is 
d2V/dv2  +(a-2qcos2v)V  =  0  (2) 

and  the  radial  (modified)  Mathieu’s  differential  equation  is 
d2U  i  du2  -  ( a~2qcosh2u)U  =  0  (3) 

where  <?  =  (&f)2  =qR+jql9  k  =o)y[ju^  and  ex  =  e\  -  jex . 

It  can  be  shown  that  there  is  a  countable  set  of  characteristic 
values  of  a,  for  which  the  solutions  are  n  or  2n  periodic. 
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One  can  easily  distinguish  between  four  kinds  of 
characteristic  values  (eigenvalues)  [15]: 
a  =  a2r(q) ,  even  solution,  period  % 

a  =  a2r+]  0?)  >  even  solution,  period  2n 
a  =  ^2r+i  (*?)  j  odd  solution,  period  n 
a  ~  ^2r+\  (4)  *  odd  solution,  period  2tt 
r  =  0,1,2,... 

For  every  characteristic  value,  the  solution  of  (2)  can  be 
given  as  a  Fourier  series: 

Se2r  (v,  q)  =  £  4m  cos  2otv  (4a) 


5e2r+1  (v,  <?)  =  £  4mt'i  cos(2/w  +  l)v 

m=0 

5o2r+2(v,^)  =  |;4:;22sin(2m  +  2)v 


So2r+l  (v,  9)  =  Z4m+i  sin(2m  +  l)v 


(4b) 

(4c) 

(4d) 


Where  the  ^’s  and  s  are  complex  Fourier  coefficients.  For 
computing  the  Fourier  coefficients,  one  has  to  obtain  first 
the  characteristic  values. 

Substitution  of  the  series  (4a),  (4b),  (4c),  and  (4d)  in  the 
differential  equation  (2),  yields  four  sets  of  equations  for 
computing  the  Fourier  coefficients: 


a2r  -q  0 
-  2q  av  -  4  -q 
0 


0  a2r  ~(2m)2 


4r 


=  o 


(5a) 


<*2,+i  -<7 

-q  ®2r+\  ~~ 


0 

-q 

-q  a2r+x-{2m  +  \f 


A™ 

A;r+' 


(5b) 


2-16 


-q  b2r+2~{2m  +  2)2 


B]r +2 


A2::; 


(5c) 


Cl-1-?  -9  0 

“ 

'£t2r+1‘ 

-9  *2,.i  -9  -9 

0 

0 

0 

-0  kr+t“(2m  +  2)2_ 

n2r+l 
.  2m+1 . 

(5d) 


We  shall  now  show  the  computation  of  a2r c  in  (5b).  All 
other  characteristic  values  can  be  calculated  in  exactly  the 
same  way.  It  is  clear  from  (5b)  that  a2r-i  can  be  seen  as  the 
eigenvalues  of  the  infinite  tridiagonal  matrix, 


I  +  q  q  0 
q  9  q 
0 


0 

q 


0  q  (2/n  +  1) 


(6) 


In  calculating  the  eigenvalues  and  Fourier  coefficients,  it  is 
necessary  to  truncate  the  matrix  after  a  finite  number  of 
columns  and  rows  to  obtain  the  required  accuracy. 

Once  the  coefficients  A£,  Jane 

determined  using  the  routine  DEVCCG  from  Microsoft 
IMSL  Math  Libraries,  the  angular  Mathieu  functions  in 
equations  (4)  as  well  as  their  respective  derivatives  can  be 
calculated.  Next,  we  compute  the  radial  Bessel  functions  of 
the  first  and  second  types  using  equation  (3)  along  with  the 
previously  computed  Fourier  coefficients.  The  equations 
used  for  the  computation  of  the  even  and  odd  functions  of 
the  first  type  are 

Jt2r(u,q)  =  (-\yJFl2  z(-l )mA22rmJm(X])Jm(x2)/A20r  (7) 


M  _  , 

X(-l)m  ^2m+l  IVm+1  (*i  Vm  (xi )  +4  (*i  )C  (*2 )]  /  4r+l 

0 

jolr(u,<i)=ny^2 

Z  )j„m]  /  B\r 

m- 1 

■/.>*,  («.?)  =  (-  i)rVw2 

Znr  4:::  tc>  <*,  y.  c )  -  jm  y  W+1  (*2 » /  4 


(8) 

(9) 

(10) 


Jq  Jq 

where  xx^—e\x2  =  —e~\  By  replacing  Jm  by  Ym, 

we  obtain  the  even  and  odd  radial  functions  of  the  second 
type  as  follows 

yc2r («. <?) = (- or s(-ir 4;l (*, y. (*2 ) /  4r  on 


Y'2r+](u,q)  =  (-VJFl2 

M 


i  (-i)m  c:  tc,  (*,  y.  (*2 )  +4  (x,  ym+1  c*2 )]  /  4 


rBir(«,q)  =  (-iYfFJ2 

t  (-i)m  4ac,  (x,  y,_,  (*2)  -  4.,  (x,)n,.,(x2)]  /  4 


A/ 

£  (-  r  4;:i  rc.  (*,  y.  (*2 )  -y.  (*,  y„,  (*2  w  /  4' 


(12) 


(13) 


(14) 


Numerical  results  for  the  derivatives,  with  respect  to  <7,  of 
the  even  and  odd  angular  and  radial  Mathieu  functions  are 
obtained  using  the  equations  given  in  [14-17]. 
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3.  Electromagnetic  Scattering  by  a  Lossy  Dielectric 
Elliptic  Cylinder 

Consider  the  case  of  a  linearly  polarized  electromagnetic 
plane  wave  incident  on  a  lossy  dielectric  elliptic  cylinder  (of 
permeability  p  ,  permitivity  8  and  conductivity  a  )  at  an 
angle  fa  measured  with  respect  to  the  X  axis  as  shown  in 
[25].  The  electric  field  component  of  the  TM  polarized 
plane  wave  of  amplitude  E0  is  given  by 

Eiz=E0eJk°pco**-*')  (15) 

where  k0  is  the  wave  number  in  free  space.  The  incident 

electric  field  may  be  expressed  in  terms  of  Mathieu 
functions  as  follows 

K  =  lAXl(c^)sem{C<i,i1)+ 

^  (16) 

m= 1 

JZZ 

Aem  =  E0r^^Sem(c0, cos^)  (17) 

om  Jy  em  (c0  )  om 


Nem{c)=  J [Sem(c,Jj)]2dv  (18) 

om  q  om 

where  £  =  cosh  u  and  t]  =  cos  v ,  c0  =k0F  ,  F  is  the 
semifocal  length  of  the  elliptical  cross  section,  Sem , 
while  Som  are  the  even  and  odd  angular  Mathieu  functions  of 

order  m9  respectively,  and  R^  are  the  even  and  odd 
radial  Mathieu  functions  of  the  first  kind,  and  Nem  and 
Nom  are  the  even  and  odd  normalized  functions. 

Electric  field  expansion  inside  and  outside  the  lossy 
dielectric  cylinder  can  be  obtained  by  solving  the  Helmholtz 
equation  in  the  elliptical  coordinate  system  using  the 
separation  of  variables  technique.  The  scattered  electric 
field  outside  the  elliptic  cylinder  for  (£>£l)  can  be 

expressed  in  terms  of  a  sum  series  of  Mathieu  and  modified 
Mathieu  functions  as  follows 

K  =  f.Ben,RlAJ(c0,OSem(c0,J?)  + 

”=0  (19) 

iBml£)(c0,dS'm(c0,ti) 

m-\ 

where  Bem  and  Bom  are  the  unknown  scattered  field 
expansion  coefficients,  R{e„  and  R ^  are  the  even  and  odd 
modified  Mathieu  functions  of  the  fourth  kind.  Similarly, 
the  transmitted  electric  field  inside  the  cylinder  for  £  < 
can  be  written  as 


K  =  I  C„R™(cx,S)S„(cvrj)  + 

m=0  (20) 

Z  ComR^(c^)Som(c„T}) 

m= 1 

where  cx  = kxF  is  complex,  and k}  =  G)^pix£x  . 

=€l-j€i  ,  Cernmd  C0OTare  the  unknown  transmitted 
field  expansion  coefficients.  The  magnetic  field  component 
inside  and  outside  the  cylinder  can  be  obtained  using 
Maxwell’s  equations,  i.e. 

». -^T^-  (21) 

6)jUxh  8v 

—  /  dE 

H „  =-J-r^r  (22) 

0)jU}h  du 

h  =  F^ cosh2  u  -  cos2  v  (23) 

The  unknown  expansion  coefficients  can  be  obtained  by 
imposing  the  continuity  of  the  tangential  field  components 
at  surface  of  the  cylinder  £  =  ,  i.e. 


I  Me„Jcuca)Bem 


Cj£uj) 


=  Z  Menm(c^)Aem  ~ ) 

w=0  Ren 


X  M«,m(Cl.C0)S„ra  firRW  (c0,£) 


=  Z  Mo„m(C],C0)4OT  oA) 

(C0 » C1 )  =  Rm  (C, ,  l)Se„  (C0 , 7)tfr 


where  fir  is  the  relative  permeability  of  the  dielectric 
region.  The  TE  scattering  by  a  lossy  elliptic  cylinder  is  also 
presented.  In  this  case  a  derivation  dual  to  that  of  the  TM 
case  leads  to  the  following  solution  [25] 

I  MenJc„c0)B™  erR%\c0,&)- = 

\>4i) 


X  Menm(c\,c0)Aem  j£>(c 
”=°  Re,,(c i,5) 


C(C,^1) 


I  Monm(cvc0)BZ 


Ron 


MAom 


m=°  l  Rj(c  nft)  j 

The  scattered  near  and  far  fields  for  the  TM  and  TE  cases 
can  be  calculated  once  the  scattered  fields  expansion 
coefficients  are  known.  The  far  scattered  field  expressions 
can  be  obtained  as  follows 
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e:  = 


J-e-*»'Ljm  [B™Sem(c0,Tfi  +  B™Som(c0,n)](  29) 
HP 


H 


"  =  Jr-e~MP  S  Jm  fcx  (co  ,  7)  +  BlEmSom (c0 , 77)]  (30) 
Po/7 

Far  Field  data  are  usually  expressed  in  terms  of  the 
scattering  cross  section  per  unit  length,  Le.,  the  echo  width. 
For  the  TM  polarization  case  it  is  defined  as 

\e:\2 

® tm  =  27tp  lim  -  —  (31) 

|£-| 

Eq.  (30)  can  be  put  in  simpler  form  for  the  purpose  of 
computation  as  follows 

A  -  =  X r  [BemScm (c0 ,  rf)  +  BomSnm  (c0 ,  ?)]  (32) 

\  A  m 


4.  Numerical  Results 

To  determine  the  accuracy  of  the  computer  program,  the 
Fourier  expansion  coefficients  of  Mathieu  functions  are 
computed.  These  results  are  in  excellent  agreement  with 
those  given  in  [16]  for  the  case  of  real  positive  argument  q , 
and  the  maximum  difference  was  found  to  be  10*5.  Also,  to 
verify  the  accuracy  of  the  computer  program  for  the  case  of 
complex  argument  q ,  we  solved  the  problem  of 
electromagnetic  scattering  by  a  lossy  dielectric  elliptic 
cylinder  since  no  results  are  available  in  the  literature.  The 
computed  numerical  results  were  compared  with  the  special 
case  of  scattering  by  a  lossy  circular  dielectric  cylinder,  and 
with  electromagnetic  scattering  by  weakly  lossy  multilayer 
elliptic  cylinder  [23,26].  Fig.  1  shows  the  numerical  results 
of  the  scattered  normalized  echo  width  a  of  a  lossy 
dielectric  circular  cylinder  having  a  relative  permittivity  of 
£r  —\  —  yl  1 .3  and  electrical  radius k0 a  =  3.33  for  TM  and 
TE  polarizations.  These  results  compare  very  well  with 
those  presented  in  [26].  In  this  case  the  circular  cylinder 
was  approximated  using  an  axial  ratio  close  to  1,  and  was 
excited  by  an  incident  plane  wave  with  <f>t  =  1 80° . 

Fig.  2  compares  the  back  scattering  normalized  echo  width 
k0<j(0)!4  for  a  homogeneous  elliptical  cylinder  with  axial 

ratio  alb-  2  and  £  =4.0.  Two  cases  are  studied;  these 
are  the  lossless  cylinder  case  with  £  =  0.0  and  a  weak  loss 

cylinder  case  where  £  =  0.6  .  The  calculated  results  are  in 
very  good  agreement  with  those  published  in  [23].  One  can 
easily  recognize  the  reduction  in  the  back  scattering 
normalized  echo  width  due  to  the  presence  of  a  small  loss. 
This  can  be  attributed  to  the  fact  that  the  weak  loss  material 
absorbed  part  of  the  incident  wave  energy.  Another  clear 
observation  is  the  shift  in  of  resonance  location.  For  a  free 
loss  material,  resonance  is  observed  around  k0a  =3.6,  while 
in  the  weak  loss  case  the  resonance  has  shifted  to  around 


k0a=2A.  Fig.  3  illustrates  the  results  obtained  in  [16]  for 
the  radial  Mathieu  function  of  the  first  kind  and  zero  order 
with  respect  to  u  for  different  values  of  q ,  and  presented 
here  to  validate  the  accuracy  of  the  computer  code.  Finally, 
Figs.  4-13  present  a  collection  of  data  for  Mathieu  functions 
(radial  and  angular)  of  complex  arguments.  The  Fourier 
expansion  coefficients  of  Mathieu  functions  for  different 
arguments  are  tabulated  in  Tables  1-4.  These  results  should 
be  very  helpful  for  researchers  solving  problems  in  the 
elliptic  coordinate  system. 

5.  Conclusions 

The  theory  for  calculating  the  angular  and  radial  Mathieu 
functions  of  complex  arguments  has  been  outlined.  The 
computed  backscattered  echo  width  of  a  lossy  dielectric 
elliptic  cylinder  was  compared  with  the  special  case  of  lossy 
circular  and  weakly  lossy  elliptic  cylinders  and  the  results 
are  in  complete  agreement.  Selected  numerical  results  for 
Mathieu  functions  are  plotted  and  tabulated  for  limited 
ranges  due  to  limitation  on  space. 
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Figure  1:  TM  and  TE  Normalized  backscattered  echo  width 
versus  <p  of  a  lossy  circular  cylinder  with  k0a  =  3.33 and 

er=l-j\l.3. 
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Figure  2:  TM  Normalized  backscattered  echo  width  versus 
<j)  of  a  lossy  dieletric  elliptic  cylinder  with  er  =  4 ,  a/b=2 

and  <j)i  =  0° 


Figure  3:  Radial  Mathieu  function  of  the  first  kind  and  zero 
order  vs  u  for  different  values  of  q. 


Figure  5:  Imaginary  part  of  radial  Mathieu  function  of  the 
first  kind  and  zero  order  vs  u  for  different  values  of  q. 


Figure  6:  Real  part  of  the  radial  Mathieu  function  of  the  first 
kind  and  zero  order  vs  u  for  different  values  of  q . 


Figure  7:  Imaginary  part  of  the  radial  Mathieu  function  of 

Figure  4:  Real  part  of  the  radial  Mathieu  function  of  the  first  the  first  kind  and  zero  order  vs  u  for  different  values  of 

kind  and  zero  order  vs  u  for  different  values  of  complex  q.  complex  q. 
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Figure  8:  Real  part  of  the  radial  Mathieu  function  of  the  first 
kind  and  zero  order  vs  u  for  different  values  of  q. 


Figure  9:  Imaginary  part  of  the  radial  Mathieu  function  of 
the  first  kind  and  zero  order  vs  u  for  different  values  of  q. 


Figure  1 1 :  Imaginary  part  of  even  periodic  Mathieu  function 
with  orders  0-5  and  q  =  5  -  j2  vs  v. 


Figure  12:  Real  part  of  odd  periodic  Mathieu  function  with 
orders  1-5  and  q  =  5  -  j5  vs  v. 


1.4 


Figure  10:  Real  part  of  even  periodic  Mathieu  function  with 
orders  0-5  and  q  =  5- j2  vs  v. 


Figure  13:  Imaginary  part  of  odd  periodic  Mathieu  function 
with  orders  1-5  and  q  =  5  -  y'5  vs  v. 
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Table  1 :  Complex  Fourier  expansion  coefficients  of  even- 
even  Mathieu  function  (  ) 


m\r 

0 

2 

g  =  5  +  jO 

Real 

Imag. 

Real 

Imag. 

0 

0.540612446 

0 

0.438737166 

0 

2 

0 

0.65364026 

0 

4 

0.14792709 

0 

-0.426578935 

0 

6 

0 

0.075885673 

0 

q  =  5  +  j  5 

0 

0.49810151 

0 

0.314442843 

0.171487341 

2 

•0.662750285 

-0.130176082 

0.604942355 

0 

4 

0.168770876 

0.133382581 

0.026159729 

-0.59020582 

6 

ns 

-0.140969223 

0.089859201 

<7  =  104-  y‘10 

0 

-0.44951258 

0.056260383 

0.029961136 

2 

0.694591697 

0 

-0.059497099 

0.255064311 

4 

-0.296630326 

0.746844887 

0 

6 

0.058531831 

0.064924524 

BBS 

Table  3:  Complex  Fourier  expansion  coefficients  of  odd- 
even  Mathieu  function  ( B Ym ) 


m\r 

2 

10 

q  =  5+j0 

Real 

Imag. 

Real 

Imag. 

2 

0.933429442 

0 

0.000033444 

0 

4 

-0.354803915 

0 

0.000642976 

0 

6 

0.052963729 

0 

0.010784806 

0 

8 

-0.004295885 

0 

0.13767512 

0 

q  =  5  +  j5 

2 

0.870294564 

0 

-0.000130329 

m 

-0,396630117 

-0.272694975 

-0.001280155 

0.001222013 

6 

0.028211436 

0.099835145 

BM 

0.021019319 

8 

0.0053059 

0.132565552 

0.136516812 

<7=10+ yio 

2 

0.734255554 

0 

-0.001900098 

-0.000156527 

-0.580164625 

-0.25100546 

mi 

0.008264692 

6 

0.137051442 

0.19916337 

0.076604876 

8 

0.002015162 

-0.052591626 

0.231256474 

0.260037365 

Table  2:  Complex  Fourier  expansion  coefficients  of  even- 
odd  Mathieu  function  ( l, ) 

j2r+\ 


m\r 

1 

3 

q=5+J0 

Real 

Imag. 

Real 

Imag. 

1 

0.762463687 

0 

0.077685798 

0 

3 

-0.63159632 

0 

0.30375103 

0 

5 

0.139684806 

0 

0.927728396 

0 

7 

-0.014915596 

0 

-0.201706148 

0 

q  =  5  +  j5 

1 

-0.517275706 

0.286011205 

-0.041658559 

0.131675807 

3 

0.765785524 

0 

0.239854654 

0.290697818 

5 

■0.197812838 

-0.153451695 

0.878673205 

0 

7 

0.007517743 

0.03807304 

-0.166961766 

■0.194677341 

<7=io+ yio 

1 

-0.433393417 

0.130050198 

0.666818493 

0 

3 

0.775016392 

0 

0.494260893 

-0.299695541 

5 

-0.373034841 

-0.201569455 

-0.01872188 

-0.426395143 

7 

0.054220182 

0.106990887 

-0.184358865 

0.054869005 

Table  4:  Complex  Fourier  expansion  coefficients  of  odd-odd 
Mathieu  function  (  ) 


m\r 

1 

3 

q  =  5  +  j0 

Real 

Imag. 

Real 

Imag. 

l 

0.940019022 

0 

0.050382462 

0 

m 

mum 

0 

0.297365513 

0 

5 

HUB 

0 

0.931566996 

0 

H 

0 

0 

<7  =  5  +  y5 

| 

1 

0.904373155 

0 

0.016113831 

0.092448918 

3 

RRm 

9RH1 

0.268912553 

0.288134908 

5 

0.057923929 

0.074976514 

0.87655541 

0 

a 

BBS 

-0.16648902 

-0.196783024 

<7  =  10  +  ylO 

a 

0.852110863 

0 

0.203330605 

3 

-0.460992769 

-0.18202088 

0.234907973 

0.510680316 

5 

0.083408677 

0.142179121 

0.741126923 

0 

7 

0.009774806 

-0.031241645 

iilH 

-0.285889059 
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Abstract:  Comprehensive  shielding  in  modern  electronic 
equipment  may  lead  to  resonant  behaviour  within  the 
equipment  enclosure.  This  paper  presents  a  method  for 
optimising  the  placement  of  sources  of  electromagnetic 
(EM)  energy  and  susceptors  to  this  EM  energy  within  an 
enclosed  resonant  cavity.  The  source  and  susceptor  are 
placed  on  a  dielectric  slab  within  a  perfectly  conducting 
enclosure  to  reduce  the  level  of  EM  coupling  between  the 
two.  Optimisation  is  facilitated  using  a  genetic  algorithm 
coupled  wTith  a  finite-difference  time-domain  solver.  Re¬ 
sults  are  presented  for  single  objective  optimisation  and 
multi-objective  optimisation  cases. 

1  Introduction 

The  explosive  increase  in  the  use  of  electronic  equipment 
in  the  information  age  has  led  to  electromagnetic  com¬ 
patibility  (EMC)  and  electromagnetic  interference  (EMI) 
issues  becoming  more  important  to  designers  of  elec¬ 
tronic  equipment.  These  issues  must  be  considered  at  the 
time  of  design  and  not  once  designs  are  at  the  prototype 
stage.  Higher  clock  speeds  and  faster  switching  transi¬ 
tions  lead  to  greater  levels  of  electromagnetic  emissions, 
while  higher  component  integration  and  lower  power  de¬ 
mands  lead  to  greater  sensitivity.  So  at  once  systems  are 
becoming  more  prone  to  generating  and  also  more  sensi¬ 
tive  towards  EMI. 

An  area  of  interest  to  the  authors  is  electronic  enclo¬ 
sures  and  the  placement  of  devices  inside  these  enclosures. 
Modern  electronic  items  have  many  different  sources  of 
EMI.  These  sources  are  often  placed  inside  a  rectangu¬ 
lar  shaped  metal  box  to  limit  the  EM  emissions  from  the 
equipment.  Care  is  taken  to  ensure  the  shielding  is  as 
comprehensive  as  possible  thus  apertures  are  kept  to  a 
minimum  and  covered  in  metallic  blanking  plates.  Such 
a  shielded,  rectangular  enclosure  has  a  good  chance  of 
forming  a  resonant  cavity  so  that  field  strengths  gener¬ 
ated  by  the  circuit  may  well  be  enhanced  once  the  cir¬ 
cuit  is  mounted  within  the  casing.  Normally  the  circuit 


may  have  a  number  of  sources  within  the  enclosure,  all 
of  which  add  to  the  EM  fields  within;  changing  the  po¬ 
sitions  of  the  various  sources  will  change  the  amount  of 
resonance  and  constructive  interference,  [1].  There  are 
also  likely  to  be  a  number  of  susceptors  at  multiple  loca¬ 
tions  throughout  the  circuit. 

There  is  obviously  an  optimal  size  and  shape  of  enclo¬ 
sure  and  component  layout  however,  exhaustively  placing 
components  inside  enclosures  and  then  computing  the  re¬ 
sultant  fields  is  a  task  that  that  would  require  a  massive 
undertaking  on  behalf  of  the  designer.  A  far  better  ap¬ 
proach  to  component  placement  is  to  use  some  kind  of 
optimisation  method  to  place  the  components  to  achieve 
a  desired  radiation  level.  This  paper  describes  the  novel 
use  of  a  genetic  algorithm  (GA)  and  a  finite-difference 
time-domain  (FDTD)  solver  to  optimise  source  and  sus¬ 
ceptor  placement  inside  a  perfectly  conducting  structure, 
building  upon  initial  work  completed  in  [2] . 

Genetic  algorithms  are  briefly  introduced  in  section  2  and 
then  the  Finite-Difference  Time- Domain  method  is  intro¬ 
duced  in  section  3.  Section  4  describes  the  setup  of  the 
computer  simulation  and  section  5  discusses  the  results  of 
these  simulations.  Finally,  section  6  presents  conclusions 
from  the  work. 

2  Genetic  Algorithms 

GAs  are  a  stochastic  search  mechanism  with  their  oper¬ 
ation  firmly  rooted  in  natural  selection  and  survival  of 
the  fittest,  [3]  and  [4].  GAs  have  been  shown  to  pro¬ 
vide  robust  search  and  optimisation  in  complex  spaces, 
[3]  and  [5].  They  use  simple  operations  on  a  population 
of  individuals,  which  lead  to  an  emergent  evolution  of  an 
individual  or  individuals.  Each  individual  in  the  popula¬ 
tion  represents  a  potential  solution  to  the  given  problem 
scenario  and  as  such  is  evaluated.  After  an  individual 
has  been  evaluated  a  figure  of  merit  (FoM)  is  attributed 
to  the  individual.  This  FoM  is  a  measure  of  how  well  the 
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Figure  1:  GA  Process  Flowchart. 

individual  solves  the  problem.  GAs  can  lead  to  the  opti¬ 
mal  solution  to  a  problem,  more  often  however,  they  are 
used  to  optimise  a  problem  towards  an  improved  solution 
facilitating  a  trade-off  between  excessive  computational 
time  and  meaningful  results.  The  GA  process  can  be 
summarised  graphically  as  shown  in  figure  1. 

The  GA  methodology  used  in  this  work  is  the  micro  ge¬ 
netic  algorithm  (jj, GA).  The  fiG A  technique  has  been 
shown  to  reach  the  optimal  area  of  multi-modal  solution 
spaces  earlier  than  conventional  GA  methods,  with  min¬ 
imal  computing  time,  [6].  It  has  already  been  demon¬ 
strated  that  the  fiGA  can  be  successfully  linked  with 
an  FDTD  solver,  [7],  for  optimisation  in  EM  problems. 
Initially  in  a  GA  a  fixed  size  population  is  created  and 
populated  with  randomly  generated  individuals  of  a  fixed 
length,  this  length  is  determined  by  the  encoding  of  the 
problem  parameters.  The  population  size  for  the  /iGA  is 
generally  kept  small  e.g.  five.  This  is  different  from  the 
classic  GA  where  the  population  size  is  typically  much 
larger,  in  the  order  of  100s.  These  individuals  are  all 
possible  solutions  to  the  same  given  problem.  They  are 
assessed  and  each  one  given  a  score,  the  FoM,  which  is 
then  used  to  determine  the  likelihood  of  reproduction  into 
the  next  generation.  The  FoM  is  analogous  to  fitness  in 
the  natural  world  where  fitter  individuals  survive  longer 
and  hence  have  a  better  chance  of  continuing  their  ge¬ 
netic  lineage.  The  particularly  successful,  higher  scoring 
individuals  may  be  reproduced  more  than  once  into  the 
next  generation.  In  the  /xGA  used  here  tournament  se¬ 
lection  is  employed.  This  is  perhaps  the  most  effective 
for  many  application  types  as  it  has  been  shown  to  pro¬ 
vide  better  convergence  towards  a  solution  in  the  initial 
stages  of  optimisation,  [8].  Tournament  selection  works 
by  randomly  choosing  N  members  of  the  population,  and 
in  this  instance,  as  in  many  others,  N  =  2.  These  indi¬ 
viduals  are  then  pitched  against  each  other  to  determine 
which  has  the  better  FoM,  the  winner  of  the  tournament 
being  selected  to  be  a  member  of  the  new  population;  this 
is  repeated  until  the  new  population  is  filled.  In  conjunc¬ 


tion  with  tournament  selection  elitist  reproduction  is  also 
employed.  This  guarantees  that  the  best  individual  from 
a  population  is  present  in  the  next  population,  hence  pre¬ 
serving  the  best  current  solution  to  the  given  problem  and 
maintaining  a  good  genetic  stock.  After  reproduction  the 
individuals  undergo  crossover.  Here  two  randomly  picked 
individuals  are  mated  together,  swapping  information  be¬ 
tween  their  chromosomes.  In  a  classic  GA  mutation  is 
also  applied,  however  in  the  /xGA  mutation  is  turned  off. 
Mutation  serves  to  alter  the  genetic  makeup  of  these  new 
offspring  with  a  small  fixed  probability.  Once  again  mim¬ 
icking  nature,  mutation  can  lead  to  either  a  detrimental 
or  beneficial  effect  on  performance.  Mutation  allows,  in 
a  classic  GA,  random  search  to  take  place  and  hopefully 
leads  to  the  avoidance  of  local  optima.  Once  the  popu¬ 
lation  has  converged  to  a  determined  optimal  value,  fur¬ 
ther  GA  operating  is  no  longer  needed  and  the  GA  can 
be  halted.  Due  to  the  small  population  size  used  in  the 
pGA  premature  convergence  is  often  seen  to  be  happen¬ 
ing.  To  prevent  this  a  restart  mechanism  is  used.  This 
restart  mechanism  involves  checking  for  convergence  in 
the  current  generation  of  individuals  and  then  restarting 
the  next  generation  with  only  the  elite  individual  and  new 
randomly  created  individuals.  The  use  of  a  restart  oper¬ 
ator  ensures  random  search  takes  place  and  leads  fiG A 
away  from  local  optima.  Convergence  is  checked  for  in 
this  application  by  measuring  the  changes  in  the  chromo¬ 
somes,  the  individuals,  between  generations.  Once  there 
are  few  changes  between  generations  then  the  population 
can  be  considered  to  be  converged;  the  setting  of  this 
limit  is  a  parameter  that  the  user  has  control  over.  An 
important  component  of  any  GA  code  is  the  device  by 
which  random  numbers  are  generated,  or  to  be  more  cor¬ 
rect  pseudo-random  numbers.  This  is  usually  a  portable 
pseudo-random  number  generator  (PRNG)  that  produces 
the  same  sequence  of  random  numbers  for  a  given  seed 
value.  The  quality  of  the  PRNG  is  an  important  factor. 
The  area  of  PRNG  research  is  vast  and  not  going  to  be, 
for  brevity,  discussed  here.  Common  references  on  the 
subject  are  [9],  [10],  [11]  and  [12]. 

A  clear  difference  can  be  seen  here  between  optimisation 
based  on  calculus  methods  and  GA  based  optimisation. 
GAs  make  use  of  interim  performance  in  the  optimisation 
problem;  calculus  based  methods  are  only  concerned  with 
optimisation  toward  an  overall  optimal  point  and  do  not 
“remember”  any  interim  performances. 

3  The  FDTD  Method 

The  FDTD  method,  [13],  [14],  has  become  one  of  the  most 
popular  methods  for  solving  Maxwell’s  equations.  It  is  a 
volumetric  domain  decomposition  technique  that  is  sec¬ 
ond  order  accurate  in  space  and  time,  easily  implemented 
in  software  and  accurately  models  the  physical  world.  It 
is  a  widely  used  method  for  EMC/ EMI  work,  [15],  as 
well  as  radar,  bioengineering  and  antenna  analysis.  The 
FDTD  method  is  described  widely  in  the  literature,  [16], 
and  so  only  a  very  brief  description  is  given  here. 

Maxwell’s  equations,  in  differential  form,  equations  (1)  - 
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(2),  are  simply  modified  to  central-difference  equations, 
discretised,  and  implemented  in  software.  The  equations 
are  solved  in  a  leap-frog  manner,  [14] ,  i.  e.  the  electric  field 
is  solved  at  a  given  instant  in  time,  then  the  magnetic  field 
is  solved  at  the  next  instant  in  time,  and  the  process  is 
repeated  over  and  over  again.  In  equations  (1)  and  (2) 
H  and  E  have  their  usual  meanings. 


9H 

9t 


(-V  x  E  +  -h), 
\V  V  ) 


(1) 


dE 

dt 


-V  x  H  -  -E. 
e  e 


(2) 


Using  a  three  dimensional  Cartesian  co-ordinate  system 
we  can  now  write  out  the  vector  components  of  the  curl 
operator.  Given  below  is  only  one  of  these  components, 
namely  the  x  component  of  the  H  field, 
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where  (i  is  magnetic  permeability  and  pf  is  a  magnetic  loss 
parameter.  It  is  the  complete  set  of  coupled  partial  dif¬ 
ferential  equations,  six  in  total,  that  are  the  fundamental 
basis  of  the  numerical  FDTD  algorithm.  The  discretised 
form  of  equation  (3)  is  shown  below: 
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where  n  is  the  time  step  under  consideration  and  i,j,  k 
are  the  three  orthogonal  spatial  co-ordinates.  The  integer 
array  MEDIA  defines  material  conditions  for  each  field 
vector  component.  This  allows  the  medium  at  each  point 
to  be  mapped  out.  The  arrays  Da,Db  are  magnetic  field 
update  coefficients,  there  are  corresponding  electric  field 
update  coefficients  also.  A  full  explanation  of  all  of  these 
equations  is  presented  in  [13]. 


4  Simulation  Setup 

The  initial  results  are  obtained  from  using  the  GA  to 
place  a  source  and  susceptor  point  relative  to  each  other 
on  the  surface  of  a  dielectric  slab  (DS)  to  achieve  the  low¬ 
est  electromagnetic  coupling  between  the  two.  The  com¬ 
plete  simulation  is  implemented  in  Fortran  77  compiled 
on  a  SUN  ENTERPRISE  server  with  512MB  of  RAM 
utilising  4  of  8  processors. 


Figure  2:  Enclosure  geometry  showing  the  representations 
of  the  internal  structure  and  dielectric  slab ,  all  sizes  are 
in  millimetres,  mm. 

4.1  Physical  Problem  Description 

The  problem  geometry  is  that  of  a  simplified  metal  box, 
which  has  perfectly  electrically  conducting  (PEC)  walls, 
an  internal  structure  and  DS  representation.  The  inter¬ 
nal  structure  is  modelled  as  two  PEC  sheets  in  the  top 
corner  of  the  geometry.  The  DS  can  be  thought  of  as  a 
representation  of  the  substrate  of  a  printed  circuit  board 
(PCB).  The  problem  geometry  is  shown  in  figure  2.  The 
DS  has  a  relative  permittivity,  tr,  of  4  and  possesses  unity 
magnetic  permeability  pr. 

4.2  n GA  and  FDTD  Setup 

The  pGA  used  is  based  on  the  implementation  by  Carroll, 
[17],  modified  to  accommodate  the  given  task  of  moving 
the  source  point  and  susceptor  point  on  the  DS.  The  two 
points  are  specified  in  a  three  dimensional  Cartesian  sys¬ 
tem  with  the  y  ordinate  being  held  constant  as  this  rep¬ 
resents  the  surface  of  the  DS.  Binary  encoding  is  used  in 
the  pGA  and  this  leads  to  a  chromosome  string  length 
of  24  bits,  4  bits  per  x,  y  and  z  co-ordinate  of  the  source 
and  susceptor.  It  should  be  noted  that  it  would  be  pos¬ 
sible  to  omit  the  y  ordinate  from  the  chromosome;  how¬ 
ever,  the  memory  saving  would  be  insignificant  especially 
when  compared  to  the  coding  effort  required  to  convert 
the  GA  for  this  one  specific  case.  Care  must  be  taken 
to  avoid  the  y  ordinate  being  altered  during  reproduction 
and  crossover,  this  is  achieved  using  a  uniform  crossover 
operator.  Uniform  crossover,  [18],  has  also  been  found  to 
give  faster  convergence  than  single  point  crossover  in  a 
pGA,  [6],  [19].  The  population  size  is  maintained  at  five. 
These  co-ordinate  values  are  passed  to  the  FDTD  solver 
which  computes  the  resulting  field  distribution  inside  the 
enclosure.  The  peak  electric  field  at  the  susceptor  point 
is  returned  as  the  FoM  to  the  pGA,  no  fitness  scaling  is 
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used  as  tournament  selection  is  the  method  of  selection 
employed  in  this  GA  implementation,  [18].  The  //GA 
attempts  to  minimise  this  FoM,  i.e.  minimise  the  peak 
electric  field  at  the  susceptor  point. 

The  FDTD  solver  is  based  on  the  equations  given  in  [13]. 
The  DS  is  one  cell  thick  and  has  one  cell  of  free  space  be¬ 
tween  itself  and  the  enclosure  wall.  The  PEC  sheets  are 
modelled  as  infinitely  thin.  A  soft  Ez  directed  sinusoidal, 
continuous  wave,  ideal  Hertzian  dipole  of  2.05 GHz  fre¬ 
quency  and  amplitude  of  10Vm-1  is  used  as  the  source 
in  the  FDTD  simulation.  The  domain  is  discretised  into 
dimensions  of  7.2727mm  in  the  x  and  2  directions  and 
7.2mm  in  the  y  direction.  Based  on  these  dimensions 
and  the  material  within  the  computational  domain  the 
time  step  size  for  the  solver  is  set  at  the  Courant  limit, 
in  this  instance  approximately  14 ps. 

Initially  two  optimisation  cases  are  run;  a  non-lossy  and 
then  a  lossy  dielectric  slab,  in  the  lossy  case  the  conduc¬ 
tivity  is  set  to  0.045m-1.  After  this  a  multi-objective 
simulation  is  run  where  two  objectives  are  to  be  op¬ 
timised,  these  being  the  reduction  in  coupling  between 
source  and  susceptor  points  at  two  specific  frequencies. 

5  Results 

5.1  Single  Objective  Optimisation 

The  progress  of  the  optimisation  process  can  be  seen  in 
figure  3  where  the  results  of  both  the  non-lossy  and  lossy 
simulations  are  presented.  This  figure  shows  the  evolu¬ 
tion  of  the  best  individual  in  the  simulation,  the  elite 
individual.  While  this  individual  remains  the  best  the 
FoM  remains  at  the  same  value,  when  a  superior  one  is 
bred,  a  step  change  occurs  indicating  the  presence  of  a 
new  elite  value. 


Figure  3:  Generation  optimisation  of  both  the  non-lossy 
and  lossy  cases . 

The  final  value  for  the  non-lossy  case  is  0.0541Vm-1 


from  an  initial  value  of  0.1172Fm-1,  and  the  final 
value  for  the  lossy  case  is  0.0240 Vm~l  from  a  value  of 
O.OfiSSFm-1.  Both  of  the  curves  are  normalised  to  the 
value  of  0.1172 Vm-1  for  ease  of  comparison.  As  expected 
the  values  of  field  strength  for  the  lossy  case  are  lower 
than  those  for  the  non-lossy  case  as  in  the  lossy  case  en¬ 
ergy  is  dissipated  in  the  dielectric  slab.  Figure  3  also 
shows  us  that  the  lossy  case  reaches  its  best  FoM  at  gen¬ 
eration  66,  earlier  than  the  non-lossy  case  which  reaches 
its  best  FoM  at  generation  91.  This  is  a  random  dynamic 
of  the  GA  and  attributable  to  the  solution  surface  topol¬ 
ogy  and  the  choice  of  PRNG  used  with  the  //GA.  The 
number  of  generations  was  limited  to  100  as  experience 
with  this  code  has  shown  that  there  is  a  minimal  return 
on  computing  time  by  exceeding  this  point.  Testimony 
to  this  is  given  by  the  fact  that  running  the  code  to  200 
generations  gives  only  a  final  value  of  0.0507V^m-1  for 
the  non-lossy  case  and  0.0204Vrm-1  for  the  lossy  case,  a 
marginal  increase  in  performance  for  a  doubling  of  com¬ 
putation  time.  This  embraces  one  philosophy  of  using  a 
GA,  namely  to  produce  an  acceptable  improvement  to¬ 
ward  an  optimal  point  for  a  limited  amount  of  effort. 

The  final  field  distribution  on  the  surface  of  the  DS  can 
be  seen  in  figure  4  for  the  non-lossy  case  and  in  figure  5 
for  the  lossy  case.  These  plots  are  generated  by  applying 
a  peak  hold  to  the  Ez  component  at  each  location  on  the 
surface  of  the  DS  throughout  the  simulation. 

The  final  positions  of  the  source  and  susceptor  points  on 
the  DS  axe  indicated;  for  the  non-lossy  case  these  are  at 
(2,16)  and  (3,35)  respectively,  and  for  the  lossy  case  the 
co-ordinates  are  (2,35)  and  (23,11).  These  co-ordinate 
values  given  are  of  course  on  the  x  and  z  axes  as  the 
y  ordinate  is  constant.  These  false  colour  plots  clearly 
show  that  the  //GA  has  indeed  found  good  solutions  but 
not  the  best  solutions.  The  global  minimum  point  in  fig¬ 
ure  4  is  at  (2,2)  with  a  value  that  is  2.1  dB  down  on  the 
position  found  by  the  //GA.  In  figure  5  the  global  mini¬ 
mum  is  at  (12,14)  with  a  value  that  is  lAdB  down  the 
//GA  position.  The  final  values  found  by  the  //GA  are 
57 .6dB  and  6A.6dB  down  on  the  global  maximum  for  the 
non-lossy  and  lossy  cases  respectively.  It  should  be  noted 
that  the  patterns  presented  in  these  figures  are  ones  of 
over  1  million  possible  patterns  due  to  source  and  sus¬ 
ceptor  placement.  Finding  the  global  minimum  on  each 
of  them  would  require  direct  evaluation  of  each  pattern. 
The  //GA  vastly  cuts  down  on  the  number  of  evaluations 
required  to  arrive  at  an  acceptable  result. 

5.2  Multi-objective  Simulations 

Multi-objective  optimisation  (MO)  is  often  sought  in  prac¬ 
tice  as  often  compliance  in  one  particular  objective  upsets 
the  compliance  of  another  objective.  MO  aims  to  achieve 
a  solution  that  can  not  be  improved  upon,  simultaneously, 
in  each  of  the  objectives.  This  is  called  a  Pareto  optimal 
solution,  and  the  set  of  all  these  solutions  is  called  the 
Pareto  optimal  set,  [20].  For  the  MO  simulations  some 
changes  had  to  be  made  to  the  GA-FDTD  code.  As  two 
objectives  were  being  optimised  a  method  of  measuring 
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Figure  4:  False  colour  plot  of  final  field  distribution  on 
surface  of  dielectric  slab  after  non-lossy  simulation.  The 
x  and  z  axes  are  marked  in  FDTD  cells .  Source  and 
susceptor  positions  are  indicated  on  the  figure. 

their  fitness  at  one  and  the  same  time  is  required.  To  ac¬ 
complish  this  the  method  of  objective  weighting  is  used, 
[20] .  This  method  is  explained  mathematically  by  equa¬ 
tion  (5), 

N 

Z  =  ^T^Wifi(x)  where  x  £  X.  (5) 

t=i 

In  equation  (5)  Z  is  a  scalar  valued  variable  that  is  a 
weighted  sum  of  the  individual  objectives,  where  there 
are  N  objectives  in  total.  The  function  /  is  the  function 
that  returns  the  FoM,  in  this  case  the  FDTD  solver,  and 
x  represents  the  parameters  of  the  function,  co-ordinate 
values  in  this  case.  The  feasible  region  of  solutions  is 
denoted  by  X.  The  sum  of  the  individual  weights  adds  to 
1  and  each  weight  lies  in  the  range  0  to  1  i.e.  0  <  Wi  <  1. 
It  is  the  scalar  value  Z  that  is  optimised  by  the  GA.  The 
method  of  objective  weighting  is  an  easily  implementable 
scheme  that  produces  a  solution  from  the  Pareto  optimal 
set. 

The  problem  setup  for  the  MO  optimisation  involves  the 
same  geometry  as  previously  with  the  only  change  being 
the  source.  A  Gaussian  pulse  modulated  onto  a  sine  wave 
is  now  used  as  the  source.  This  gives  a  symmetrical  spec¬ 
trum  about  the  carrier  frequency,  no  DC  component,  and 
a  bandwidth  controlled  by  the  length  of  the  pulse.  Its 
mathematical  form  given  in  [13]  is  described  by  equation 
(6), 

r  "-"p  ia 

Ez  =  Eoe  L "Jecov  J  sin(2irfo(n  -  no) At).  (6) 

The  Gaussian  pulse  is  centred  around  frequency  fo  at 
step  no  with  a  1/e  characteristic  decay  of  ndecay  steps, 


Figure  5:  False  colour  plot  of  final  field  distribution  on 
surface  of  dielectric  slab  after  lossy  simulation.  The  x  and 
z  axes  are  marked  in  FDTD  cells.  Source  and  susceptor 
positions  are  indicated  on  the  figure . 

and  At  is  the  step  size.  The  time  series  and  frequency 
response  of  the  source  is  shown  in  figure  6.  The  centre 
frequency  for  the  source  is  chosen  at  1.500 GHz  and  the 
objectives  for  optimisation  were  chosen  as  the  minimisa¬ 
tion  of  the  1.002 GHz  and  2.004 GHz  components  in  the 
frequency  response  which  is  computed  on  the  fly  as  the 
FDTD  simulation  progresses  as  detailed  in  [13].  These 
frequencies  were  chosen  as  they  tie  in  with  equipment 
used  in  a  related  measurement  study.  Only  non-lossy  sim¬ 
ulations  were  run  in  this  setup. 

The  resulting  frequency  response  from  the  Ez  component 
at  the  susceptor  point  the  simulations  can  be  seen  in  fig¬ 
ure  7.  The  dashed  line  in  figure  7  shows  the  initial  DFT 
before  the  optimisation  process  begins;  the  solid  line  in 
the  figure  displays  the  final  DFT  after  completion  of  the 
optimisation  process.  A  clear  difference  in  the  two  curves 
at  the  respective  objective  frequencies  can  be  seen.  Also, 
it  can  be  clearly  seen  that  the  field  values  at  the  objective 
frequencies,  on  the  final  DFT,  are  considerably  lower  than 
at  other  frequencies,  within  the  same  response,  excepting 
of  course  at  the  tails  of  the  response  where  there  is  mini¬ 
mal  energy  from  the  source.  The  values  of  the  objectives 
on  the  final  response  are  54 dB  down  for  the  1.002 GHz 
objective  on  the  max  value  in  the  response  and  the  value 
at  2.004 GHz  is  44 dB  down.  Relative  to  the  initial  DFT 
the  final  response  in  down  27 dB  at  the  1.002 GHz  objec¬ 
tive  and  down  35 dB  at  the  2.004 GHz  objective. 

6  Conclusions 

A  method  of  using  a  /iGA  in  conjunction  with  a  FDTD 
solver  to  facilitate  electromagnetic  optimisation  has  been 
shown.  The  conjoining  of  the  two  codes  allows  a  dif¬ 
ficult  design  problem  to  be  tackled,  namely  that  of  ra- 
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Frequency  (Hz) 


Figure  6:  Time  series  of  the  pulsed  source,  upper  graph , 
and  its  DFT,  lower  graph .  This  is  the  source  character¬ 
istic  used  in  the  multi-objective  optimisation  simulations. 

diating  and  susceptible  component  placement  within  an 
arbitrary  metallic  structure  to  minimise  coupling.  The 
pGA  technique  is  used  as  this  gives  a  significant  reduc¬ 
tion  in  computation  time  with  little  loss  in  performance 
of  the  final  optimised  result.  Ensuring  that  the  pGA  finds 
the  global  optimum  in  difficult  optimisation  problems  is 
an  area  that  needs  further  attention,  but  the  ability  of 
the  technique  to  reach  a  near  optimal  solution  has  been 
demonstrated.  Both  single  and  multiple  objectives  for 
optimisation  have  been  presented  with  promising  results 
from  each.  It  is  worth  noting  that  if  an  exhaustive  search 
were  to  be  completed  on  the  problems  presented  then  over 
one  million  simulations  would  be  required  taking  over  one 
year  to  complete;  the  GA  based  search,  however,  took  ap¬ 
proximately  four  hours. 

The  design  optimisation  cases  presented  are  rather  sim¬ 
plistic  but  do  prove  the  concept  of  the  hybrid  code.  More 
challenging  problem  geometries  can  easily  be  accommo¬ 
dated  in  the  FDTD  code  with  little  burden  to  comput¬ 
ing  time  as  once  the  domain  has  been  discretised  it  is  the 
number  of  resulting  cells  that  chiefly  determines  the  com¬ 
putation  time.  It  is  envisaged  that  this  tool  can  be  used 
to  provide  design  rules  for  component  placement  within 
an  enclosed  resonant  cavity  and  that  it  can  can  also  be 
used  to  assess  specific  designs. 
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ABSTRA  CT.  The  aim  of  this  paper  is  to  present  an  approach , 
able  to  deal  with  all  possible  connections  of  voltage  and 
current  sources  and  impedances ,  combining  conductors  in 
which  the  skin  effect  is  taken  into  account  and  conductors  in 
which  skin  effect  is  neglected This  approach  is  obtained 
using  a  bottom-up  methodology .  In  this  way,  the  meaning  of 
terms  in  the  generalized  approach  is  naturally  inherited  from 
some  specific  problems .  This  model  is  presented  in  a  compact 
form ,  preserving  sparse,  symmetric  and  positive-definiteness 
matrices .  The  vectors  and  matrices,  computed  during  the 
solution  stage ,  are  employed  in  engineering  calculations  as 
current  and  inductance  computations  providing  compact 
expressions  suitable  for  efficient  algorithms.  Finally,  the 
proposed  approach,  implemented  in  a  FEM  software  package 
developed  by  the  authors,  is  applied  to  the  study  of  a  three- 
phase  transformer,  supplied  with  a  balanced  three-phase 
voltage  (sinusoidal  and  nonsinusoidal)  and  loaded  with  an 
unbalanced  three-phase  RL  impedance.  The  agreement 
between  the  computed  and  experimental  results  shows  the 
validity  of  the  proposed  model  and  its  implementation. 

1  INTRODUCTION 

The  natural  excitation  for  the  Finite  Element  Method  (FEM) 
in  electromagnetism  is  current  (current  density).  This  current 
is  supplied  by  an  external  electric  circuit;  and,  usually,  the 
value  of  this  current  is  unknown,  as  the  circuits  consist  of 
passive  elements  and  voltage  sources.  Therefore,  FEM  and 
electrical  equations  are  coupled  by  these  unknown  currents. 
There  are  two  different  approaches  to  couple  these  equations. 
In  the  first  approach,  called  direct  coupling,  the  equations 
obtained  from  FEM  and  the  electrical  circuits  are  integrated  in 
a  single  system  of  equations  having  the  magnetic  vector 
potential  as  its  solution.  In  the  second  one,  called  indirect 
coupling,  two  different  systems  of  equations  are  obtained;  one 
of  them  takes  into  account  the  FEM  equations,  and  the  other 
one,  the  circuit  equations.  In  this  paper,  a  method,  based  on 
the  direct  coupling  of  FEM  and  circuit  equations  is  presented. 

Many  approaches  on  coupled  2D  Finite  Element  models  and 
circuit  equations  have  been  proposed  f  l]-[7].  Most  of  them, 
take  into  account  the  skin  effect  for  all  the  conductors  or,  they 
neglect  its  influence,  without  evaluating  its  magnitude  over 
each  conductor.  In  this  paper,  a  bottom-up  methodology  is 


applied,  generalizing  the  meaning  of  the  studied  examples, 
instead  of  a  theoretically  obtained  approach.  So,  some 
common  winding  connections  were  studied—one-phase  circuit 
and  three-phase  circuits  with  star  and  delta  connections—, 
being  presented  here  the  three-phase  star  connection  handling 
unbalanced  load  and  voltage  conditions.  In  this  paper,  the 
model  is  applied  to  the  study  of  a  power  transformer,  although 
it  has  been  applied  to  other  non-static  machines  as  induction 
motors,  [10]. 


2  BASIC  EQUATIONS 

Assuming  that  the  displacement  currents  are  neglected, 
Maxwell’s  equations  using  a  magnetic  vector  potential 
formulation  can  be  written  in  2D  as: 


=  -  J(x,  y)  (1) 


Where  A(x,y)  is  the  magnetic  vector  potential,  J(x,y)  the 
current  density  (z  component)  and  v  the  reluctivity  that 
depends  on  the  magnetic  flux  density  vector  (B).  In 
anisotropic  materials,  reluctivity  depends  on  |B|  and  the  angle 
(0)  between  B  and  the  rolling  direction  of  the  material 
v=v(zi,e). 

The  expression  of  current  density  depends  on  how  the 
conductors  are  defined.  So,  the  conductors  can  be  modelled  as 
solid  conductors,  in  which  case  the  skin  effect  is  taken  into 
account,  or  stranded  conductors,  in  which  case  the  skin  effect 
is  neglected.  Stranded  conductors  are  particularly  useful  to 
model  windings  in  which  an  individual  mesh  considering  their 
relations  would  exceed  in  many  cases  the  requirements  of 
actual  computers  without  providing  a  significantly  greater 
accuracy.  If  all  the  conductors  in  the  problem  are  solid 
conductors,  the  following  equation  is  employed,  [8]  [9]: 

r)\ 

V x  vV  x  A  =  <rEeit  -  a—  -crvc  x  V  x  A  (2) 


The  first  term  on  the  right  side  in  (2),  E ^  being  the  external 
electric  field  and  a  the  electrical  conductivity,  represents  the 
current  density  due  to  the  voltage  source.  The  second  term 
takes  into  account  the  current  density  due  to  the  temporal 
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variation  of  A.  The  third  one  represents  the  current  density 
due  to  the  movement,  being  ve  the  velocity  of  the  conductor 
with  respect  to  B.  By  using  a  frame  of  reference,  such  that  the 
relative  velocity  becomes  zero,  (2)  can  be  written  as,  [8]: 

dA 

V  x  vV  x  A  =  crEexf  -  cr —  (3) 

dt 


The  total  current  of  a  conductor  can  be  obtained  by  integrating 
the  right  side  of  (3): 


In  2D  problems  (3)  becomes: 

-V(\A7A)-ar  —  -<j—  (5) 

/  dt 

Where  us  is  the  voltage  between  the  terminals  of  the  conductor 
and  /  the  length  of  the  conductor  (thickness  of  the  considered 
domain  along  the  z-axis).  Equation  (5)  takes  into  account  the 
skin  effect. 

When  the  conductors  are  considered  stranded,  it  is  supposed 
that  the  current  (if)  is  the  same  in  all  the  conductors  connected 
in  series  (winding);  then,  their  current  density  can  be  obtained 
by  J^Nif/S.  Where  S  is  the  entire  section  of  all  the  conductors 
(N turns)  integrating  the  winding,  a  the  electrical  conductivity 
and  Rf  the  DC  resistance.  The  current  can  be  calculated  in  2D 
problems,  by  considering  the  voltage  K/r  between  the  terminals 
of  the  winding  as  the  sum  of  the  voltages  of  its  conductors: 


Using  (1)  and  (6): 


Figure  1.  The  windings  are  star  connected. 

connected  (figure  1).  They  are  modelled  with  FEM,  and 
indeed  can  or  can  not  share  the  same  2D  magnetic  circuit.  The 
total  DC  resistances,  including  the  resistances  corresponding 
to  the  windings  and  the  load,  are  called  Rr,  Rs,  Rj.  In  the 
same  way  LR,  Ls,  LT  take  into  account  the  inductance 
corresponding  to  the  load  and  the  part  of  leakage  flux  not 
considered  in  the  two-dimensional  model. 

From  figure  1  the  following  equations  can  be  obtained: 

v*  =  +  (*,'■,  +  4  J]  -  [Rrh  +  4^~]  ( 1  °) 

=eT-eR + (/Mr + 4  ^)  -  [y, +  4^f)  o  » 

iR  +  is  +  iT=0  (12) 


(7) 


Summarizing  (3)  and  (7)  into  a  single  equation  for  2D 
problems: 


3  EXAMPLE.  THE  WINDINGS  ARE  STAR 
CONNECTED 

The  authors  have  solved  some  complex  winding  connections. 
In  this  case,  three  windings  of  stranded  conductors  are  star 
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Being  (w}T  the  phase  vector,  which  provides  information 
about  the  topology  of  each  winding,  (p(x,y)  the  shape  or 
interpolation  functions,  whose  formulas  are  well-known  for 
classical  elements,  and  Q  the  entire  domain  where  the  integral 
is  extended  by  using  the  phase  functions  ft/. 


Pj(x,y) 


l 

-l 

o 


If  (. x9y )  €  s; 
If  (x,y)  €  5; 

Otherwise 


(15) 


Using  the  Crank-Nicholson  time  stepping  method,  the 
following  equations  are  obtained: 


Using  (20),  (21)  and  the  definition  of  Am  in  (14)  the  following 
equation  is  obtained: 

((X-H4)  <22> 

Where  the  term  [  W]  is  called  the  phase  matrix  defined  as  the 
matrix  containing  the  phase  vectors  presented  in  the  problem 
(three  in  this  case). 

Using  (1)  and  (22): 


(16)  ([5,+i']  +  /[lF][^][Gi][^][lFf){^}  = 

(17)  -MM[Ga][G,]{/} 


Using  (9),  (10),  (1 1),  (12),  (13),  (16)  and  (17): 


Ms+Mt  ~Mt 
Mt  Mt  +  Mr 
Ms  Mr 


l 

Mi 


where  [St+At]  is  the  stiffness  matrix.  Therefore,  equation  (23) 
models  2D  magnetic  circuits  with  a  three  star  connected 
windings  of  stranded  conductors,  the  magnetic  vector 
potential  in  the  nodes  at  t+At  being  the  unknown  values. 
Besides,  current  in  the  windings  at  t  must  be  computed  using 
efficient  algorithms. 

4  THE  GENERALIZED  SYSTEM  OF 
EQUATIONS 


M'r(Ms  +  Mt) 

M  'sMt 

~M  'tMs 

1  1 

M'rMt 

M's(Mt+Mr) 

M'tMr 

1  is  1 

The  example  analyzed  above  can  be  generalized  to  different 

M  'rMs 

M'sMr 

M1t(Ms  +  Mr)_ 

{ \ 

types  of  connections  and  conductors  (solid  and  stranded)  [11]. 

where: 

Ms  =  RjtSt  +  2Lj  j  =  R,SJ 

M'j  =  RjAt-2Lj  j  =  R,S,T  (19) 

M\  =  MrMs  +  MSMT  +  MtMr 

Equation  (18)  can  be  written  in  a  packed  form  as: 


(24> 

-a^iwo.iio.i-ww-i 


[g.1M(R"}-{X})-[g.][g,]{,'} 

The  current  density  can  be  obtained  using: 


0 

0 

’s;*' 

Sr 

Vs 

c 

,/^A/ 

lR 

Js* 

►  = 

0 

0 

is 

sr 

0 

* s 

0 

■s^ 

i _ 

.t+At 

lT 

So,  by  example,  if  a  specific  problem  includes  windings  of 
(20)  stranded  conductors  with  star  connection  and  other  windings 
with  whatever  connections,  the  contribution  of  the  star 
connected  windings  to  the  terms  of  (24)  is  obtained  from  (23). 

The  matrices  in  (24)  maintain  the  sparsity,  symmetry  and 
positive-definiteness  of  the  traditional  FEM.  Obviously, 
sparsity  depends  on  geometry  and  the  number  of  nodes  since 
the  mean  magnetic  vector  potential  introduces  connections 
(2U  among  the  nodes  of  the  windings  modelled  with  stranded 
conductors  (14). 

Defining  «,  m  and  p  as  the  number  of  nodes,  and  windings 
made  of  stranded  and  solid  conductors  respectively,  the 
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meaning  and  dimensions  of  the  terms  in  (24)  can  be  expressed 
as: 


•  [S]nxn  is  the  stiffness  matrix. 

•  [7]nxn  is  called  the  mass  matrix,  containing  information 
about  the  current  distribution. 

•  is  the  phase  matrix  showing  information  about  the 
distribution  of  windings  made  of  stranded  conductors.  It  is 
composed  of  the  phase  vectors  defined  in  (14). 

•  [Ga]mxm,  [Gb]mxm  and  [Gc]mxm  take  into  account  how  the 
windings  defined  in  [W]  are  connected.  In  (18)  and  (20) 
they  are  defined  to  star  connections. 

•  is  the  vector  of  unknown  nodal  values  representing 
the  value  of  the  magnetic  vector  potential  at  the  current 
time  step  t+At 

•  {A*}n  represents  the  value  of  the  magnetic  vector  potential 
at  the  previous  time  step  t. 

•  [Q]  has  the  same  meaning  in  solid  conductors  as  [W]  in 
stranded  conductors. 

•  {unbelt1  \p  rePresents  the  addition  of  voltages  applied  to 

the  solid  conductors  at  the  current  time  step  t+At  and  at  the 
previous  time  step  t. 

•  \u{^)+t  represents  the  addition  of  voltages  applied  to 

the  stranded  conductors  at  the  current  time  step  t+At  and  at 
the  previous  time  step  t. 

•  \?cab$m  represents  the  value  of  current  at  the  previous  time 
step  in  conductors  defined  as  stranded. 

5  COMPUTATION  OF  ENGINEERING 
RESULTS 

As  shown  in  (23)  and  (24),  currents  in  windings  must  be 
calculated  during  the  solution  stage.  Besides,  this 
computational  cost  must  be  kept  as  low  as  possible.  Other 
engineering  results  such  as  inductance,  magnetic  flux  and 
electromotive  force  can  be  obtained  by  making  use  of  the 
terms  defined  in  (24),  computed  at  each  time  step  avoiding 
time  consuming  calculations  in  the  postprocessing  stage. 


Current  in  windings  made  of  conductors  considered  as 
stranded  ones  can  be  calculated  using  (1),  (5)  and  the 
definition  of  phase  vector: 


[5'+A']{^+A'}  =  {/'+A'}  =  X^{^}  if*  (27) 

k-1 


Where  m  is  the  number  of  windings.  The  current  of  a  winding 
(0  can  be  computed  multiplying  both  terms  in  (27)  by  the 
transposed  i-phase  vector  {a>i)T,  taking  into  account  that  the 
product  of  {o>i}T  and  {a>k}  is  null  when  k  is  not  the  same  as  i: 


n,  {a,}7  {<*>,} 


(28) 


This  method  computes  the  phase  current  using  as  data,  the 
phase  vector,  the  number  of  turns,  together  with  the  stiffness 
matrix  and  the  vector  of  the  magnetic  vector  potential  nodal 
values  at  the  current  time  step  t+At.  Therefore,  all  the  terms  in 
(28)  are  computed  during  the  solution  stage  of  FEM. 


5.2  Flux  and  Electromotive  Force  Calculations 

These  results  can  be  easily  computed  using  (13)  and  (14). 


5.3  Inductance  Calculations 


Inductance  has  a  geometrical  meaning  when  all  the  materials 
exhibit  a  linear  behaviour.  In  other  case,  several  non¬ 
coincident  definitions  can  be  taken  into  account  [12].  Using 
the  definition  based  on  flux  linkages,  self-inductance  can  be 
calculated  in  windings  made  of  stranded  conductors  as: 


L'+a'  = 


N±£*_  IN'toflA1'*} 


.  /-*•  A/ 

U 


(29) 


5.1  Current  Calculations 

The  winding  currents  can  be  calculated  using: 

i/t)=  jj(t)  dQ 

D/ 


In  solid  conductors  (25)  can  be  expressed  as: 


In  the  same  way  mutual  inductance  can  be  computed  using  the 
flux  linkage  definition: 


(25) 


,  M  t  +  &t  —  Nl 
IvJ  ij - 


.  1+ A/ 

h 


{g£[g^l  {a™) 


Nj  {<y,f  {coj} 


(30) 


(26)  As  the  current  is  computed  in  (28),  self  and  mutual 
inductances  are  calculated  using  the  phase  vectors  and  the 
stiffness  matrix,  that  can  be  dependent  on  time.  Stiffness 
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matrix  is  affected  by  the  B-H  curve  when  it  exhibits  a 
nonlinear  behaviour. 

In  [13]  the  inductance  computation  in  case  of  solid  conductors 
is  stated. 

6  EXAMPLES 

A  software  package  has  been  developed  implementing  (24).  It 
has  been  tested  with  many  examples  of  increasing  complexity 
over  transformers,  considering  several  balanced  and 
unbalanced  connections  in  primary  and  secondary  -one  phase 
and  three  phase  transformers-,  steady  and  unsteady  states, 
sinusoidal  and  nonsinusoidal  voltages  and  current  sources 
(balanced  and  unbalanced);  two  of  these  examples  are 
presented  below. 

6.1  Steady-state 

The  model  is  applied  to  a  3  kVA  three-phase  AY  transformer 
(figure  2).  The  transformer  is  supplied  with  a  balanced 


Load 

Figure  2.  The  transformer  is  delta-star  connected. 


nonsinusoidal  three-phase  voltage  (figure  3),  generated  with  a 
programmable  AC  power  source.  These  waveforms  are 
applied  to  the  software  package  using  files  obtained  from  an 
oscilloscope.  The  secondary  is  loaded  with  an  unbalanced 
three-phase  RL  impedance.  Figure  4  shows  the  simulated 
currents  in  steady  state,  and  figure  5  shows  the  computed  and 
measured  currents. 


0  0.005  0.01  0.015  0.02  0.025  0.03  0,035  0.04 


0  0.005  0.01  0.015  0.02  0.025  0.03  0.035  0.04 


0  0.005  0.01  0.015  0.02  0.025  0.03  0.035  0.04 


Secondary  phase  currents  -Time  (s) 

Figure  4.  Primary  and  secondary  computed  currents. 


00°5  0.01  0.015  0.02  0.025  0.03  0  0.01  0.02  0  0.01  0.02  0  0.01 

Time  (s)  Time  (s)  Tlme(s)  Time  (s) 


Figure  3.  Line-to-line  supplied  voltages. 


Figure  5.  Computed  (solid)  and  measured  (dashed)  currents. 
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LAB/Primary  line  currents  -  Time  (s) 

Figure  6.  Computed  (FEM)  and  measured  (LAB)  currents. 


6.2  Non  steady-state 

The  transformer  is  now  YY  connected,  and  loaded  with  the 
same  unbalanced  three-phase  RL  impedance.  Figure  6  shows 
the  computed  and  measured  primary  currents;  the 
discrepancies  between  them  are  due  to  the  unfavourable 
conditions  of  the  high  supply  voltage  and  the  instant  chosen 
for  the  switching-on  of  the  transformer. 

7  CONCLUSIONS 

A  generalized  approach  for  FEM  and  external  circuits 
coupling  has  been  presented  considering  both  types  of 
conductors  simultaneously,  stranded  and  solid  respectively. 
This  model  has  been  established  using  a  bottom-up 
methodology.  The  meanings  of  the  terms  of  the  particularized 
equation  (23)  have  been  generalized  to  obtain  the  terms  in 
general  equation  (24).  Indeed,  when  a  problem  containing  a 
particularized  example  is  selected,  the  expressions  of  these  are 
included  in  the  terms  of  (24). 

Current,  magnetic  flux,  electromotive  force  and  inductance 
definitions  have  been  obtained  using  some  terms  of  (24), 
simplifying  their  notation  and  speeding-up  their  computation 
during  the  solution  stage  rather  than  during  the  postprocessing 
stage. 

Finally,  the  approach,  implemented  in  a  software  package 
developed  by  the  authors,  has  been  applied  to  several  cases. 
The  results  have  shown  a  good  agreement  with  those 
measured  in  the  lab. 
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Abstract-  Characteristic  modes  developed  by 
Garbacz,  Harrington  and  Mautz  have  long  been 
used  in  the  analysis  of  radiation  and  scattering 
from  conducting  bodies  and  apertures.  For  their 
computation,  it  is  necessary  to  solve  an  eigen- 
system  of  the  form  Ax  —  ABx.  If  the  matrices 
(A,B)  are  Hermitian  and  B  is  positive  definite, 
the  generalized  eigenvalue  problem  can  be  accu¬ 
rately  solved  using  the  simultaneous  diagonal¬ 
ization  technique  (SDT).  Because  of  numerical 
approximations  and  rounding  sometimes  it  may 
happen  that  the  matrices  properties  deteriorate 
and  the  SDT  procedure  becomes  inapplicable. 
In  this  work  a  new  technique,  developed  recently 
by  Higham  and  Cheng  is  proposed  as  a  method 
to  solve  these  deteriorate  cases.  It  is  applied  to 
the  computation  of  the  characteristic  modes  for 
some  scattering  problems.  Results  are  analyzed 
and  discussed. 


1.  Introduction 

The  generalized  eigenvalue  problem: 

Ax  =  ABx  (1) 

is  often  encountered  in  computational  electro¬ 
magnetics  [1].  The  solution  of  (1)  depends  on 
properties  of  the  matrices  A  and  B  resulting 


from  the  discretization  of  the  field  equations.  If 
these  matrices  are  not  singular,  then  (1)  can  be 
reduced  to  the  solution  of  the  standard  eigen¬ 
problem: 

Hx  =  Ax  with  H  =  B-1  A  (2) 

However,  if  A  and  B  are  Hermitian  and  B  is 
positive  definite,  the  problem  (1)  can  be  solved 
with  a  higher  degree  of  accuracy  applying  the  si¬ 
multaneous  diagonalization  technique  (SDT)  to 
the  matrices  (A,B)  than  the  direct  calculation 
of  the  inverse  matrix  B"1  [2].  Although  the 
underlined  properties  are  sometimes  held  of  the 
matrices  issued  in  electromagnetic  problems,  un¬ 
fortunately  the  discretization  performed  by  nu¬ 
merical  algorithms  on  the  electromagnetic  field 
equations  often  causes  the  loss  of  these  proper¬ 
ties.  As  a  consequence,  the  SDT  becomes  in¬ 
applicable.  Recently,  Higham  and  Cheng  have 
faced  the  above  problem  from  a  theoretical  point 
of  view  and  they  have  proposed  a  technique  that 
allows  to  use  the  SDT  even  in  presence  of  dis¬ 
cretization  errors.  This  method  is  based  on  the 
concept  of  nearest  definite  pair  [3].  In  this  work, 
we  have  applied  Higham-Cheng  results  to  evalu¬ 
ation  of  the  characteristic  modes  for  the  scatter¬ 
ing  from  conducting  bodies.  The  paper  is  orga¬ 
nized  as  follows:  section  II  presents  the  basics  on 
the  theory  of  characteristic  modes  for  conduct- 
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ing  bodies.  Section  III  and  IV  illustrate  the  SDT 
procedure  for  the  resolution  of  the  positive  defi¬ 
nite  Hermitian  generalized  eigenproblem  and  its 
extension.  In  section  V  we  show  the  numerical 
results.  Finally,  in  section  VI,  the  conclusions. 


2.  Theory  of  Characteristic  Modes  for 
Conducting  Bodies 

The  characteristic  modes  have  been  intensively 
used  in  the  analysis  of  radiation  and  scatter¬ 
ing  from  conducting  bodies  and  apertures  [4]- 
[13].  They  are  numerical  entire  basis  functions 
that  in  principle  can  be  computed  for  any  ob¬ 
ject  shape.  Since  these  eigenfunctions  include 
the  behavior  of  unknown  current  density  flow¬ 
ing  on  the  metallic  bodies  and  apertures,  only 
a  small  number  of  them  are  required  for  a  good 
reconstruction  of  it.  So,  if  available,  character¬ 
istic  modes  would  lead  to  a  scattering  or  radia¬ 
tion  problem  treatable  with  a  classical  Method 
of  Moment  even  in  presence  of  a  large  number  of 
objects  [10]-[13].  In  the  case  of  scattering  or  ra¬ 
diation  from  conducting  bodies,  these  modes  are 
basically  solutions  of  a  generalized  eigenproblem 
involving  the  impedance  operator  Z,  which  re¬ 
lates  the  surface  current  J  on  a  conducting  body 
to  the  tangential  component  of  the  incident  elec¬ 
tric  field  on  he.  [5]: 

Z3  =  4an  (3) 

One  method  to  solve  a  functional  equation  of  the 
form  (3)  is  to  obtain  a  modal  solution  in  terms 
of  eigenfunctions  of  Z  .  For  this  purpose,  we 
consider  the  operational  eigenproblem: 

ZJk  =  VkMJk  (4) 

in  which  M  is  a  suitable  weight  operator  and 
Uk  =  1  +  j\k-  The  operator  Z  can  be  uniquely 
represented  in  the  form: 


Note  that  7 Z  and  X  are  real  selfadjoint  operators 
and  furthermore  %  is  positive  definite  [5].  The 
choice  M  —  H  reduces  the  complex  operatorial 
eigenproblem  (4)  to  a  real  operatorial  eigenprob¬ 
lem: 

XJk  =  (8) 

The  eigenfunctions  Jk  have  been  named  char¬ 
acteristic  modes  [5].  They  can  be  numerically 
evaluated  through  the  reduction  of  the  opera¬ 
tor  equation  (8)  to  a  matrix  equation  using  the 
Method  of  Moment  [5], [10].  For  this  aim  the  h- 
th  modal  current  density  on  the  metallic  body 
is  expanded  in  a  set  of  N  suitable  subsectional 
basis  functions  Bn : 

JV 

j*  =  £/nSn  (9) 

n~  1 

where  the  In  are  unknown  complex  constants. 
Substitution  of  (9)  into  (8)  and  application  of 
the  Galerkin  method  leads  to  a  matrix  equation 
of  the  form: 

XJfc  =  AfcR/fc  (10) 

where  is  the  column  vector  of  the  unknown 
coefficients  In.  Equation  (10)  is  a  generalized 
eigensystem  of  the  form  (1)  in  which  the  eigen¬ 
values  A*;  and  the  eigenvectors  Ik  approximate 
the  eigenvalues  and  the  eigenfunctions  of  the 
equation  (8).  It  follows  by  X  and  7 Z  properties 
that  the  matrices  X  and  R  are  expected  Her¬ 
mitian  and  furthermore  R  is  expected  positive 
definite.  Consequently,  in  principle  it  is  possible 
to  use  the  SDT  for  solving  the  eigensystem  (1). 
But  because  of  numerical  approximation,  the  R 
matrix  often  becomes  indefinite  and  the  SDT  al¬ 
gorithm  becomes  inapplicable.  In  the  following, 
it  will  be  shown  that  the  use  of  the  nearest  def¬ 
inite  pair  technique  permits  to  employ  the  SDT 
giving  a  significant  improvement  in  the  evalua¬ 
tion  of  the  characteristic  modes  over  the  direct 
inversion  technique. 


Z  =  R  +  jX  (5) 

where  the  operators  7 Z  and  X  are  defined  as: 


3.  Resolution  of  the  Positive  Definite 
Generalized  Eigenproblem  with  the 
Simultaneous  Diagonalization  Technique 


1Z  — 

x  = 


z  +  z* 
2 

Z-Z* 
2  j 


(6)  The  generalized  eigenproblem  where  the  matri¬ 
ces  A  and  B  are  Hermitian  and  in  which  one 
is  positive  definite  plays  an  important  role  in 

(7) 

matrix  theory.  If  we  assume,  without  loss  of 
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generality,  that  B  is  positive  definite,  the  prob¬ 
lem  (1)  can  be  accurately  solved  by  means  of 
a  method  named  simultaneous  diagonalization 
technique  [2]  as  follows:  Let  Ui  be  an  unitary 
matrix  whose  columns  are  an  orthonormal  set 
of  eigenvectors  for  B.  Premultiplying  and  post- 
multiplying  (1)  by  UJ  and  Ui,  where  the  star 
denotes  the  complex  conjugate  transpose  opera¬ 
tion,  and  taking  into  account  that: 


u;ux  =  Uiu;  =  dmp(i)  (11) 

where  diag(  1)  is  the  identity  matrix,  we  obtain: 

Ay  =  AB  y  (12) 

in  which  the  matrices  A,  B  and  the  vector  y  are 
of  the  form: 

A  =  UJ  AUi 

V  =  U  \x  (13) 

B  =  UIBUi  =  diag(fi) 

and  in  which  diag(fj)  is  the  diagonal  matrix  formed 
by  eigenvalues  of  the  B  matrix.  Next,  we  intro¬ 
duce  the  nonsingular  transformation  H: 

H  ^diag(fi~2)  (14) 

Using  H  in  (12),  we  obtain: 

A  z  =  AB  z  (15) 

where  A,  z  and  B  are: 

A  a?  H*AH 

z  =  H”1^  (16) 

B  -  H*BH  -  diag(  1) 


Finally,  if  we  construct  an  unitary  matrix  U2 
whose  columns  are  an  orthonormal  set  of  eigen¬ 
vectors  for  A,  the  matrix  transformation  T  — 
U1HU2  permits  simultaneously  to  reduce  the 
matrices  A,  B  to  diagonal  form: 


T*AT  =  diag(X)  \ 
T*BT  =  diag(l)  J 


(17) 


resolving  the  eigenproblem  (1).  An  efficient  im¬ 
plementation  of  the  above  procedure  that  uti¬ 
lizes  both  the  Cholesky  factorization  and  the 
symmetric  QR  algorithm  is  described  in  [14]. 


4.  Extension  of  the  Simultaneous 
Diagonalization  Technique 

The  application  of  the  procedure  discussed  in 
the  previous  section  is  not  limited  to  the  treat¬ 
ment  of  Hermitian  pairs  (A,B)  with  B  positive 


definite,  but  sometimes  it  can  be  extended  to  the 
Hermitian  pairs  in  which  the  B  matrix  is  indef¬ 
inite.  For  this  aim,  the  concept  of  definiteness 
has  been  restated  referring  directly  to  the  ma¬ 
trix  pair  (A,  B)  rather  than  to  a  single  matrix 
B.  More  in  details,  pair  (A,  B)  having  Crawford 
number  7  [15],  denoted  as: 

7(A,B)=  min  v^A*)2  +  C?*BS)2  (18) 

M2=i 

strictly  positive  has  been  named  as  definite  and, 
when  7  =  0  the  pair  has  been  named  as  indefi¬ 
nite.  Also,  it  is  possible  evaluate  the  definiteness 
of  a  Hermitian  pair  (A,B)  drawing  its  Field  of 
Values  F( A  +  jB)  in  the  complex  plane.  The 
field  of  values  is  defined  as  the  set  of  all  the  val¬ 
ues  assumed  by  the  Rayleigh  quotients  of  the 
pair  [3]: 

F(A  +  jB)=  ?*(Azt/B)*  (19) 

The  pair  will  be  definite  if  and  only  if  zero  does 
not  lie  in  F(A  +  jB),  otherwise  it  will  be  indef¬ 
inite.  A  solution  of  (1)  in  which  B  is  indefinite, 
but  the  Crawford  number  7  of  the  pair  (A,  B)  is 
positive,  has  been  given  in  [15]  by  means  of  the 
Stewart’s  theorem  [16].  This  theorem  ensured 
the  existence  of  an  angle  0  e  [0,27r)  that  per¬ 
mits  to  define,  starting  from  a  given  pair  (A,  B) 
with  7  >  0,  a  new  pair  (A0,  B$)  in  which  B0  is 
positive  definite,  as  follows: 


A#  m  Acos 0  +  Bsinfl  1  . 

B0  —  - Asin0  +  Bcosfl  J  ^  ' 

For  this  new  pair  is  now  applicable  the  SDT. 
The  eigenvalues  A  of  the  old  pair  (A,B)  and 
the  eigenvalues  A#  of  the  new  pair  (A^B#)  are 
simply  related  by: 

A  _  Aflcosfl  -  sinfl 

Agsin#  -f  cos#  ' 

However,  it  may  happen  that  the  discretization 
errors  may  affect  the  matrices  in  a  way  that  a 
pair  (A,B)  theoretically  expected  to  be  defi¬ 
nite  becomes  indefinite  in  practice.  As  a  con¬ 
sequence  the  SDT  cannot  be  applied.  For  this 
last  case  it  can  be  useful  evaluate  the  nearest 
definite  pair  (A  -f  AA,  B  +  AB)  and  attempt 
to  resolve  the  problem  (1)  for  this  one.  Else, 
assigned  the  indefinite  Hermitian  pair  (A,B), 
we  must  find  the  nearest  Hermitian  definite  pair 
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(Ad-  A  A,  B  +  AB),  having  a  specified  Crawford 
number  7  =  8  >  0  according  to  the  distance 

[3], [17]: 

d6( A,B)  =  mm{||[AA  AB]||2  : 

7(A  +  AA,  B  +  AB)  =  (5}  (22) 

and  the  angle  0  by  means  of  which  the  nearest 
pair  can  be  simultaneously  diagonalized.  This 
problem  has  been  recently  solved  by  Higham  and 
Cheng  by  means  of  the  following  theorem  [3]: 


Figure  1:  Geometry  for  scattering  from  circular, 
elliptic  and  square  metallic  cylinders. 


Higham-Cheng’s  Theorem:  Let  A,B  €  Cnxn 
be  Hermitian  and  let  C  =  A  +  jB  and  A^  = 
A cos{<j>)  +  B sin(<j>).  Let  min(O<0<27r)  \ max  (A*) 
be  attained  at  the  angle  6  and  let  Aq  have  the 
spectral  decomposition: 

A$  =  Qdiag(nk)Q*  with  fin  <  ...  <  jjti  (23) 

If  0  G  F(C)  then: 

ds(  A,B)  =  <5  +  /xi  (24) 

If  0  g  F(C)  then: 

^(AjB)  =  max(S  +  /xi,.0)  (25) 

In  both  cases,  two  sets  of  optimal  perturbations 
in  (22)  are: 

AAi  =  Q  diag(min(-S  +  /x*,  O))Q*cos0  1 
ABi  =  Q diag(min{—6  -f  0))Q*sin^  J 

(26) 

and 

AA2  =  -~dfi(A,  3)diag(l)cos0  1  ,  , 

AB2  =  -fib(A,B)dtop(l)sinfi  J  K  } 

A  detailed  demonstration  of  this  theorem  can  be 
found  in  [3].  The  most  important  consequence 
of  the  above  results  is  that  the  SDT  can  be  gen¬ 
eralized  to  the  Hermitian  pairs  (A,  B)  having  B 
matrix  indefinite. 

5.  Numerical  Results 

In  order  to  check  the  validity  and  accuracy  of  the 
foregoing  procedure  in  the  computation  of  char¬ 
acteristic  modes  for  conducting  bodies,  numeri¬ 
cal  results  for  some  cases  are  carried  out.  Firstly, 
characteristic  modes  have  been  computed  for  the 
TMZ  scattering  from  a  perfectly  conducting  cir¬ 
cular  cylinder  and  from  an  elliptic  one,  since  for 


Figure  2:  Plot  of  the  smaller  normalized  eigenvalues 
for  the  TMZ  scattering  from  a  circular  cylinder  of 
radius  a  =  %  (Problem  Size  128  x  128), 

these  cases  characteristic  modes  are  known  in 
analytical  form  [12].  In  Fig.(l)  are  represented 
the  geometries  for  both  problems.  In  Fig. (2)  is 
reported  the  plot  of  the  smaller  eigenvalues  of 
the  R  matrix  obtained  discretizing  the  EFIE  for 
the  scattering  due  to  a  plane  wave  impinging 
on  a  metallic  circular  cylinder  with  radius 
Pulse  basis  functions  and  delta  function  test¬ 
ing  functions  are  used  to  create  the  matrix  R 
[1].  Since  only  the  properties  of  the  matrices 
above  are  considered  here  for  the  computation 
of  characteristic  modes,  it  is  not  necessary  to 
specify  the  direction  of  the  incident  field.  The 
depicted  eigenvalues  are  scaled  to  show  the  rel¬ 
ative  value  only.  This  plot  clearly  shows  that 
the  R  matrix  is  indefinite,  even  if  it  is  expected 
definite  because  the  properties  of  the  operator 
H.  In  Fig.  (3)  is  plotted  the  field  of  values  for 
the  pair  (X,R).  Because  the  origin  of  the  com¬ 
plex  plane  is  not  contained  in  F(X  +  fR)  this 
pair  is  definite.  In  Figs. (4)  and  (5)  are  reported 
the  plot  of  the  smaller  eigenvalues  of  the  R  ma- 
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Figure  3:  Field  of  Values  for  pair  related  to  compu¬ 
tation  of  characteristic  modes  for  circular  cylinder. 

trix  and  the  field  of  values  F(X  +  y’R)  for  the 
pair  obtained  considering  the  scattering  from  a 
metallic  elliptic  cylinder  with  semi  major  axis 
a  =  A  and  semi  minor  axis  b  =  ^  by  an  inci¬ 
dent  plane  wave.  For  this  scattering  problem  we 
have  obtained  both  an  indefinite  matrix  B  and 
an  indefinite  pair  (X,  R)  as  clearly  depicted  in 
these  figures.  In  Table  (1)  are  shown  the  an¬ 
alytical  eigenvalues  versus  the  numerical  ones 
for  the  previously  analyzed  scattering  problems. 
Numerical  eigenvalues  have  been  obtained  using 
both  the  Higham-Cheng  procedure  and  the  di¬ 
rect  inversion  technique.  Excellent  agreement  is 
shown  comparing  the  exact  eigenvalues  and  the 
numerical  ones  computed  by  means  of  the  former 
procedure  while  the  eigenvalues  obtained  using 
the  direct  inversion  technique  are  very  inaccu¬ 
rate.  As  last  case,  consider  the  evaluation  of 
the  characteristic  modes  for  the  scattering  from 
a  metallic  square  cylinder  of  side  a  =  A.  In 
Figs.(l)  and  (6)  are  reported  the  geometry  for 
the  problem  and  its  field  of  values,  respectively. 
As  it  is  shown  in  Fig. (6),  the  matrix  pair  (X,R) 
for  this  case  is  definite  because  the  field  of  val¬ 
ues  does  not  contain  the  origin  of  the  complex 
plane.  In  Fig.  (7)  is  shown  the  obtained  current 
distribution  on  the  cylinder  by  a  plane  wave  im¬ 
pinging  at  an  angle  <fio  =  0.  The  continuous  line 
indicates  the  current  computed  by  means  of  the 
standard  MoM  while  crosses  indicate  the  solu¬ 
tion  obtained  using  characteristic  modes  com¬ 
puted  using  the  Higham-Cheng  procedure.  A 
very  good  agreement  is  clearly  observed. 


Figure  4:  Plot  of  the  smaller  normalized  eigenvalues 
for  the  T Mz  scattering  from  an  elliptic  cylinder  with 
semi  major  axis  a  =  A  and  semi  minor  axis  b  —  - 
(Problem  Size  180  x  180). 


Figure  5:  Field  of  Values  for  pair  related  to  compu¬ 
tation  of  characteristic  modes  for  elliptic  cylinder. 


6.  Conclusions 

A  technique  developed  by  Higham  and  Cheng 
for  the  treatment  of  the  Hermitian  generalized 
eigenproblem  based  on  the  concept  of  nearest 
definite  pair  has  been  proposed  for  the  compu¬ 
tation  of  the  characteristic  modes.  Numerical  re¬ 
sults  for  some  scattering  problems  are  presented. 
In  all  the  cases,  the  outlined  procedure,  even 
if  acting  on  very  ill  conditioned  matrices,  pro¬ 
vides  very  accurate  numerical  results.  Finally, 
we  point  out  that  the  application  of  the  method 
is  not  limited  to  the  presented  cases  being  fruit¬ 
fully  applicable  to  a  wide  number  of  electromag¬ 
netic  problems  which  posses  the  discussed  prop¬ 
erties  [18]. 
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Figure  6:  Field  of  Values  for  pair  related  to  com¬ 
putation  of  characteristic  modes  for  square  cylinder 
(Problem  Size  165  x  165). 


Figure  7:  Plot  of  the  current  density  on  the  metallic 
square  cylinder.  The  continuous  line  indicate  solu¬ 
tion  computed  by  means  of  the  MoM.  Crosses  indi¬ 
cate  solution  obtained  using  characteristic  modes. 
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ABSTRACT  The  space  formed  by  the  ground  and 
ionosphere  is  known  to  act  as  a  resonator  for 
extremely  low  frequency  (ELF)  waves .  Lightning 
discharges  trigger  this  global  resonance,  which  is 
known  as  Schumann  resonances  at  the  frequencies  of 
8,  14,  21  Hz  etc .  Even  though  the  inhomogeneity  (like 
day-night  asymmetry,  local  perturbation  etc.)  is 
important  for  such  subionospheric  ELF  propagation, 
the  previous  analyses  have  been  always  made  by 
some  approximations  because  the  problem  is  too 
complicated  to  be  analyzed  by  exact  full-wave 
analysis.  This  paper  presents  the  first  application  of 
the  conventional  FDTD  method  to  such 
subionospheric  ELF  wave  propagation,  in  which  any 
kinds  of  inhomogeneities  can  be  included  in  the 
analysis,  to  be  compared  with  the  observational 
results.  We  show  the  application  of  FDTD  to  our 
problem  and  present  a  few  numerical  computational 
results  to  be  compared  with  those  by  the  pre-existing 
analysis  method. 

1.  Introduction 

Recently  there  has  been  again  a  great  interest  in 
subionospheric  ELF  waves  because  of  the  two 
important  findings.  The  first  one  is  the  suggestion  by 
Williams  [1]  that  the  global  warming  which  is  a  very 
important  issue  for  human  being,  can  be  effectively 


Fig.  1.  Configuration  of  the  Earth-ionosphere 
cavity. 


monitored  by  the  intensity  of  Schumann  resonances. 
Schumann  resonance  is  a  global  resonance 
phenomenon  in  the  Earth-ionosphere  cavity  as  shown 
in  Fig.  1,  which  is  triggered  by  lightning  discharges 
in  the  thunder-active  equatorial  region.  The 
resonance  frequencies  are  8, 14,  21  Hz  etc.  in  the  ELF 
range.  The  second  reason  is  closely  associated  with 
the  finding  of  optical  phenomena  (red  sprites  etc.)  in 
the  mesosphere  and  the  lower  ionosphere,  and  it  is 
found  that  this  mesospheric  optical  phenomenon  is  a 
source  of  strong  ELF  signals  (called  ELF  transients) 
[2].  This  ELF  transient  signal  is  one  of  the  important 
tools  for  the  study  of  those  mesospheric  optical 
phenomena  and  then  the  electrodynamic  coupling 
between  the  atmosphere  and  ionosphere. 

There  have  been  published  a  few  monographs 
dealing  with  the  subionospheric  ELF  propagation 
[3-5].  Schumann  [6]  was  the  first  to  predict  the 
presence  of  resonances  in  the  spherical  earth- 
ionosphere  cavity  and  suggested  a  mathematical 
formulation  of  the  propagation  problem  at  ELF.  A 
great  simplification  is  the  presence  of  only  a  single 
globally  propagating,  zero-order  TEM  mode  [7], 
while  all  the  higher-order  modes  attenuate  severely  at 
distances  of  several  waveguide  effective  heights. 
Despite  this  simplification,  the  complex  structure  of 
the  lower  ionosphere  imposes  an  intricate  three- 
dimensional  electrodynamical  problem  that  cannot  be 
reduced  to  practical  techniques  without  certain 
additional  simplifying  assumptions.  This  is  the 
reason  why  several  fundamental  observed  properties 
of  Schumann  resonances  cannot  be  well  explained 
[5,8],  and  please  look  at  the  forthcoming  monograph 
on  ELF  [8]  indicating  a  lot  of  unsolved  problems  in 
the  ELF  range.  The  factors  of  making  the  problem 
very  complete  are  (1)  radially  (vertically) 
inhomogeneity  of  the  ionosphere,  (2)  day-night 
asymmetry,  (3)  local  ionospheric  perturbations  etc. 

In  this  sense  it  is  highly  required  that  very 
complicated  situations  should  be  solved  in  an  exact 
way  even  by  using  the  numerical  methods.  If  the 
exact  numerical  solutions  are  available  for  any 
complicated  cavity,  it  would  be  possible  for  us  to 
examine  the  validity  of  the  previous  approximate 
solutions. 
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2.  Three-dimensional  FDTD  analysis 

The  above  discussion  has  led  to  the  requirement 
for  a  three-dimensional  finite-difference  time  domain 
(FDTD)  application  to  our  ELF  propagation  problem 
in  spherical  coodinates  with  azimuthal  dependence. 
Fortunately,  this  3D  FDTD  code  has  been  developed 
by  Holland  [9],  which  can  be  used  in  our  application. 
Fig.  2  (right  panel)  shows  the  spherical  coordinate 
system  ( r ;  6,  0)  and  location  of  a  unit  cell.  The  left 
panel  of  Fig.  2  illustrates  the  location  of  field- 
evaluation  points  on  a  unit  cell.  As  shown  in  Fig.  2, 
we  build  a  cell  entitled  I,  J  and  K  along  each  of  three 
axes  {r,  0,  and  0)  .  In  the  r  direction  I  =  1  is  the 
starting  grid  (r- a  (Earth’s  radius)  6.4Mm)  and  this 
ground  is  assumed  to  be  a  perfect  conductor.  I  =  Nr 
is  the  outermost  radial  boundary,  is  taken  to  be  r  =  a 
+  150  km  and  to  be  a  perfect  conductor.  In  the  0 
direction  J  =  1  is  the  North  Pole  {6  =  0°)  and  J  =  N0 
is  the  South  Pole  ( 0  =  180°)  .  The  coordinate  0 
indicates  the  longitude  in  such  a  way  that  K  -  1  and 
K  -  are  the  same  point,  but  the  latter  has  encircled 
the  globe.  On  this  condition  the  six  electromagnetic 
field  components  are  expressed  by  Holland  [9]. 

The  Maxwell’s  equations  to  be  solved  are 
expressed  as  follows  [10]. 


Fig.  2.  The  spherical  coordinate  system 
(r, 0,0)  and  location  of  a  unit  cell  (right 
panel).  The  left  panel  illustrates  the  location 
of  field-evaluation  points  on  the  unit  cell. 


dE 

VxH  =  £0  —  +  0E  +  Jext 
dt 


(1) 


where  E  is  the  electric  field,  H  is  the  magnetic  field, 
€0 ,  //o  are  the  dielectric  constant  and  permeability 
of  free  space,  a  is  the  conductivity  which  is  a 
function  of  space,  and  Jejr/  is  the  current  source  of  the 
system  (here  this  is  a  lightning  discharge).  Here  we 
write  down  the  field  updating  equations  in  the 
spherical  coordinates  as  follows. 


In  the  region  of  I=l~Nr_,,  J=2~Ne.]  and  K=2~N^  the  field  updating  equation  for  Er  component  is  given  by; 


£r(l,J,K)"+l/2  =  £,(UK)"-|/2  +  (c0/A /  +  a/2)2 
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sin0(j)//(t(],J,K)n  -  sin0(j  -  \)Hfi (l,  J  - 1  ,K)n 


/?(l)sin6>0(j) 

A  <p 


a  e 
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(2) 


The  corresponding  equation  for  Ee  is  as  follows. 

Ee  (1, J,  Krl/2  =  £-^~~  Ee (l,J,  K)"-|/2  +  (g0 /At  +  a/2)2 

£q /At  -f-  (7/2 

~//r(i,j,Ky,-#f(i,j>K-i)" 
sin#(./)A0 


^o(l) 

R(\)Hf  (I,J,  K)n  -  *(I  - \)H0 (I  - 1, J,  K)1 


,«  Tl 


A r 


(3) 


and  this  equation  pertains  to  I=l~Nr.i,  J=2~Ne.i  and  K=2~N^.  And  then,  the  field  updating  equation  for  E, 
component  is  as  follows. 
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^(UK)"+1/2  (I,J,K)”-1/2  +(e0/At  +  <r/2)2 

e0/At  +  (j/2 

\^)[  Ar 

//r(l,J,K)n -Hr(\,  J-1,K)"T 
Ad 

J 

for  I=2~Nr_i,  J=l~Ne_i  and  K=2~N^  Ar,  A0  and  A<j>  are  the  grid  size  in  the  r,  0,  directions,  and  At  is  time  step. 
The  field  updating  equation  for  magnetic  field  components  (Hr,  He,  H^)  are  given  as  follows. 

//,(u,Kr =^(u,k)" 

1  sing0(j  +  l)£^(l,J  +  l,K)”+l^2  — sinfl0(j)£^(l,J,K),l+^2 

/?o(l)sin0(j)  A6 

£g(l,J,K+l)w+1/2  -£g(l,J,K)”+1/2T 
A<j> 


for  I=2~N„  J=l~Ne.,  and  K=1~N!,1. 

Mo 

1  Er(l, J,K  + 1  )”+1/2  -  Er (i, J,K)"+1/2 
R{  I)  sin<?0(j)AfZ> 

J?o(l  +  l)^(Hl,J,Kr1/2-fi0(l)^(l,] 

A r 

for  I=1~Nm,  J==2~N6-1  and  K=1~NH. 

//^(i,j,k)"+1  =//,>(i,j,k)" 

/*0 

1  £(l  +  l)£g(l  +  l,J, K)”+l/2  - R{\ )Ee (I, J, K)”+1/2 

£(l)  A r 

£r(l,J  +  1,K)"+1/2  - £f (I, J,K)”+1/2  T 
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jy 


for  I=l~Nr_b  J=l-Ne.i  and  K=KN4,1. 

The  field  components  are  on  a  modified  Yee  cell  with  unit  vectors  r,  6 ,  and  and  the  corresponding 
indices  I,  J,  K.  Spatial  locations  are  given  in  term  of  positions  Ro(I),  0o(J),  <l>o(K).  Points  lying  midway  between 


these  locations  are  given  by  R(I),  0(J),  <j)(K). 

3.  Modeling  of  conductivity  profile  and 
current  source,  and  some  numerical 
examples 

In  order  to  verify  the  usefulness  of  this  3D 
FDTD  computation,  we  show  some  numerical  results 
for  a  uniform  cavity  without  any  day -night 
asymmetry,  because  we  can  compare  our  own  results 
with  the  previous  analytical  computations  [11,12]. 


However,  the  height  profile  of  the  conductivity  is 
definitely  induced  in  our  computations. 

First  we  indicate  the  height  profile  of  the 
atmospheric  conductivity,  which  is  expressed  by, 

<j(z)  =  a0exp(z!h)  (8) 

where  r  is  the  height  above  the  ground  and  h  is  the 
scale  height  ( h  =  6km).  <J0  is  assumed  to  be 
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10'14[S/m].  This  is  a  typical  profile  [10].  Next  we 
assume  a  lightning  discharge  current  in  the  form  of 
two  exponentials  by  Bruce  and  Golde  [13],  but  we 
adopt  a  model  with  much  slower  time  variation  than 
the  actual  lightning  in  order  to  validate  the  horizontal 
grid  size  of  250km. 

I(t)  =  70{exp(-70/)-exp(-100/)}  (9) 

where  t  is  in  second.  Assuming  this  current  source, 
we  made  the  FDTD  computations  in  which  the  grid 
size  in  r  is  2km,  and  the  grid  sizes  in  0  and  0are 
/r/72  (about  250km). 

Fig.  3  illustrates  the  computational  results  on  the 
temporal  evolution  of  horizontal  magnetic  field  for 
four  different  distances  (d  =  5,  10,  15  and  20  Mm, 
indicating  the  distance  in  the  6  direction)  when  the 
current  source  is  located  at  the  equator  (0=  nil,  <j) 
~  0°).  The  pulse  is  generated  at  t  =  0  and  we  can 
understand  the  propagation  of  a  pulse  in  the  v 
waveguide;  we  notice  the  direct  wave  and  its 
corresponding  antipodal  wave.  As  you  see  from  this 
figure,  with  the  increase  in  propagation  distance,  the 
higher  frequency  components  will  decay  so  that  we 
have  a  spread  in  waveform  for  larger  propagation 
distances.  The  horizontal  magnetic  field  of  the 
antipodal  wave  is  opposite  in  sign  to  that  of  the  direct 
wave,  so  that  the  signal  for  d  ~  20Mm  is  zero  because 
of  the  complete  cancellation  of  the  direct  and 
antipodal  waves.  An  important  characteristic  from 
this  figure  is  that  we  expect  severe  damping  between 
fif=5Mm  and  lOMm,  but  negligible  damping 
between  d~  lOMm  and  15Mm.  These  numerical 
computations  are  compared  with  the  corresponding 
computation  for  their  model  [11]  (the  result  is  shown 
in  Fig.  4),  in  which  there  exist  some  differences  in  the 
modeling.  The  first  one  is  the  difference  in  the  input 
signal  (source)  because  they  used  a  delta  function. 
The  sign  of  the  horizontal  magnetic  field  is  taken  to 


Fig.  4.  Computational  result  of  the  temporal 
evolution  of  horizontal  magnetic  field  for 
four  different  propagation  distances  (d=5, 
10,  15  and  20  Mm). 


Fig.  3.  Same  as  Fig.  3,  but  with  the  previous 
analytical  method[ll].  This  figure  should 
be  compared  with  Fig.  3. 

be  oppose  to  that  in  this  paper.  We  have  found  an 
extremely  good  similarity  between  these  two  figures, 
suggesting  the  usefulness  of  this  3D  FDTD  analysis 
in  subionospheric  ELF  propagation.  Another 
difference  is  definitely  due  to  the  difference  in 
propagation  constants.  They  [11]  have  used  a 
heuristic  propagation  constant  (v3  and  v2)  based  on 
the  experimental  data, 

V(/)  =  V,  -;'v2  =(/-2)/6-///100,  (10) 

where  /  is  frequency,  and  v3,  v2  are  real  and 
imaginary  part  of  propagation  constant,  respectively. 

On  the  other  hand,  we  have  used  the  atmospheric 
conductivity  profile,  Eq.(8)  as  mentioned  before.  By 
using  the  approximate  formulas  by  Greifinger  and 
Greifinger  [14],  we  can  estimate  the  real  and 
imaginary  parts  of  the  propagation  constant  for  our 
ionospheric  model.  These  results  are  shown  in  Fig.  5. 
Fig.  5(a)  refers  to  the  real  part,  while  Fig.  5(b)  refers 
to  the  imaginary  part  (or  attenuation  constant)  of  the 
propagation  constant.  The  real  parts  are  nearly 
identical  for  the  two  models,  but  we  see  a  significant 
difference  in  the  imaginary  part  (i.e.  damping  effect). 
Especially  we  notice  a  lot  of  difference  in  the 
attenuation  constant  at  higher  frequencies.  In  the 
paper  [11],  they  used  the  heuristic  propagation 
constant,  Eq.  (10),  so  that  they  have  not  indicated  any 
explicit  conductivity  height  profile.  Of  course  it  is 
theoretically  possible  for  us  to  obtain  the  conductivity 
profile  by  some  approximations  from  Eq.  (10),  but 
this  kind  of  job  is  not  so  meaningful.  In  our  next 
paper  we  will  compare  the  attenuation  from  both 
methods  under  the  same  conductivity  profile  in  order 
to  prove  the  deviation  in  Fig.  5(b)  is  caused  by 
different  computational  models  instead  of  numerical 
errors. 

Finally,  we  show  the  last  figure  for  exact 
comparison  between  ours  and  the  previous  model 
[11]  for  the  same  propagation  distance  (d=10Mm)  in 
Fig.  6. 


Frequency  [Hz] 


Frequency  [Hz] 


Fig.  5.  Comparison  of  propagation  constant,  (a)  refers  to  the  propagation  constant  (real  part),  while  (b)  the 
imaginary  part,  indicating  attenuation  constant.  Two  models  are  compared;  our  model  and 
Nickolaenko  and  Hayakawa  [11] 


4.  Conclusions  and  future  application  of 
this  FDTD  analysis  to  much  more 
complicated  models 

In  the  previous  section  we  have  shown  only  one 
example  of  our  FDTD  analysis  for  subionospheric 
ELF  propagation  and  a  comparison  with  the  previous 
analytical  solution  has  indicated  that  our  FDTD 
application  would  be  of  great  potential  in  our  future 
studies.  As  you  see  from  the  above  discussion,  the 
field  updating  mechanism  is  exactly  the  same  as  that 
of  the  high-frequency  FDTD  [16],  and  it  seems  that 
there  is  no  special  treatment  in  the  ELF  wave  analysis. 
We  can  list  what  we  will  be  able  to  do  in  this  ELF 
field  by  means  of  our  FDTD  analysis. 

As  for  the  ELF  propagation  itself  (such  as  the 
propagation  of  ELF  transients  (Q  bursts,  slow  tails)), 
we  can  study  the  reflection  mechanism  of  these  ELF 
waves  in  the  lower  ionosphere  in  details  by 
comparing  our  FDTD  results  with  the  previously 
proposed  approximations  [14,15]. 

Then,  as  for  Schumann  resonance,  there  are 


Time  (ms) 

Fig.  6.  Comparison  of  the  ELF  waveforms  for  a 
particular  d=10Mm  estimated  by  our 
method  and  previous  method  [11]. 


several  possibilities  of  the  application  of  our  FDTD 
method.  For  example,  (1)  the  full  consideration  of 
day-night  asymmetry  in  order  to  well  explain  the 
diurnal  variation  of  Schumann  resonance  intensity 
and  frequencies,  (2)  the  local  perturbation  (such  as 
the  perturbation  at  high  latitudes  (auroral  regions)). 

We  will  use  our  FDTD  analysis  for  these  above 
questions  in  order  to  improve  the  physics  coming 
from  these  subionospheric  ELF  propagation  studies. 
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